/*
 * Copyright 2014 Actelion Pharmaceuticals Ltd., Gewerbestrasse 16, CH-4123 Allschwil, Switzerland
 *
 * This file is part of DataWarrior.
 * 
 * DataWarrior is free software: you can redistribute it and/or modify it under the terms of the
 * GNU General Public License as published by the Free Software Foundation, either version 3 of
 * the License, or (at your option) any later version.
 * 
 * DataWarrior is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
 * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 * See the GNU General Public License for more details.
 * You should have received a copy of the GNU General Public License along with DataWarrior.
 * If not, see http://www.gnu.org/licenses/.
 *
 * @author Thomas Sander
 */

package com.actelion.research.chem;

import java.io.IOException;
import java.io.ObjectInputStream;
import java.io.ObjectOutputStream;
import java.io.Serializable;

import com.actelion.research.util.Angle;

/**
 * While the Molecule class covers all primary molecule information as atom and bond properties,
 * the atom connectivity and coordinates, its derived class ExtendedMolecule handles secondary,
 * i.e. calculated molecule information. Most important are the directly connected atoms and bonds
 * of every atom, information about rings and whether they are aromatic, and some atom properties
 * that depend on their near neighbors. This calculated information is cached in helper arrays and
 * stays valid as long as the molecule's primary information is not changed. High level methods,
 * e.g. getPath(), that require valid helper arrays take care of updating the cache themselves.
 * Low level methods, e.g. isAromaticAtom(), which typically are called very often, do not check
 * the validity of or update the helper arrays themselves for performance reasons. If you use low
 * level methods, then you need to make sure that the required helper array information is valid by
 * calling ensureHelperArrays().
 * Typically you will never instantiate an ExtendedMolecule, but rather use a StereoMolecule.
 */
public class ExtendedMolecule extends Molecule implements Serializable {
    static final long serialVersionUID = 0x2006CAFE;

    /**
     * To interpret a stereo center as fisher projection, all non stereo
     * bonds must be vertical and all stereo bonds must be horizontal.
     * FISCHER_PROJECTION_LIMIT is the allowed tolerance (currently 5.0 degrees).
     */
    public static final float FISCHER_PROJECTION_LIMIT = (float)Math.PI / 36;
    public static final float STEREO_ANGLE_LIMIT = (float)Math.PI / 36;   // 5 degrees

	public static final int cMaxConnAtoms = 16; // ExtendedMolecule is not restricted anymore
                                               // However, this is a suggestion for editors and other classes
    transient private int mAtoms,mBonds;
    transient private RingCollection mRingSet;

    transient private int mPi[];
    transient private int mConnAtoms[];
    transient private int mAllConnAtoms[];
    transient private int mConnAtom[][];
    transient private int mConnBond[][];
    transient private int mConnBondOrder[][];

	public ExtendedMolecule() {
		}


	public ExtendedMolecule(int maxAtoms, int maxBonds) {
		super(maxAtoms, maxBonds);
		}

	public ExtendedMolecule(Molecule mol) {
		super(mol.getMaxAtoms(), mol.getMaxBonds());
	    mol.copyMolecule(this);
		}


	/**
	 * Copies a part of this Molecule into destMol, being defined by a mask of atoms to be included.
	 * @param destMol receives the part of this Molecule
	 * @param includeAtom defines atoms to be copied; its size may be this.getAtoms() or this.getAllAtoms()
	 * @param recognizeDelocalizedBonds defines whether disconnected delocalized bonds will keep their
	 * single/double bond status or whether the query feature 'delocalized bond' will be set
	 * @param atomMap null or int[] not smaller than includeAtom.length; receives atom indices of dest molecule
	 */
    public void copyMoleculeByAtoms(ExtendedMolecule destMol, boolean[] includeAtom, boolean recognizeDelocalizedBonds, int[] atomMap) {
        if (recognizeDelocalizedBonds)
            ensureHelperArrays(cHelperRings);

        destMol.mAtomList = null;

        int atomCount = includeAtom.length;

        if (atomMap == null)
        	atomMap = new int[atomCount];

        destMol.mAllAtoms = 0;
        for (int atom=0; atom<atomCount;atom++)
            atomMap[atom] = includeAtom[atom] ? copyAtom(destMol, atom, 0, 0) : -1;

        destMol.mAllBonds = 0;
        for (int bnd=0; bnd<mAllBonds;bnd++) {
        	int atom1 = mBondAtom[0][bnd];
        	int atom2 = mBondAtom[1][bnd];
            if (atom1<atomCount && atom2<atomCount) {
            	if (includeAtom[atom1] && includeAtom[atom2])
            		copyBond(destMol, bnd, 0, 0, atomMap, recognizeDelocalizedBonds);
            	// if we keep only one atom and have a bi-polar bonds, then we need to neutralize
            	else if (mAtomCharge[atom1] != 0
            		  && mAtomCharge[atom2] != 0
            		  && (mAtomCharge[atom1] < 0
            			^ mAtomCharge[atom2] < 0)) {
            		if (includeAtom[atom1])
            			mAtomCharge[atom1] += (mAtomCharge[atom1] < 0) ? 1 : -1;
            		if (includeAtom[atom2])
            			mAtomCharge[atom2] += (mAtomCharge[atom2] < 0) ? 1 : -1;
            		}
            	}
        	}

        copyMoleculeProperties(destMol);
        destMol.mValidHelperArrays = cHelperNone;

        destMol.renumberESRGroups(cESRTypeAnd);
        destMol.renumberESRGroups(cESRTypeOr);

        if (destMol.mAllAtoms != atomCount)
            destMol.setFragment(true);

        if (recognizeDelocalizedBonds)
            new AromaticityResolver(destMol).locateDelocalizedDoubleBonds();
        }


	/**
	 * Copies a part of this Molecule into destMol, being defined by a mask of bonds to be included.
	 * Bonds, whose atoms carry opposite charges are treated in the following manner: If only one of
	 * the two bond atoms is kept, then its absolute charge will be reduced by 1.
	 * @param destMol receives the part of this Molecule
	 * @param includeBond defines bonds to be copied
	 * @param recognizeDelocalizedBonds defines whether disconnected delocalized bonds will keep their
	 * single/double bond status or whether the query feature 'delocalized bond' will be set
	 * @param atomMap null or int[] not smaller than this.getAllAtoms()
	 * @return atom map from this to destMol with not copied atom's index being -1
	 */
    public int[] copyMoleculeByBonds(ExtendedMolecule destMol, boolean[] includeBond, boolean recognizeDelocalizedBonds, int[] atomMap) {
        if (recognizeDelocalizedBonds)
            ensureHelperArrays(cHelperRings);

        destMol.mAtomList = null;

        if (atomMap == null)
        	atomMap = new int[mAllAtoms];

        destMol.mAllAtoms = 0;
        for (int atom=0; atom<mAllAtoms;atom++) {
        	atomMap[atom] = -1;
        	for (int i=0; i<mConnAtoms[atom]; i++) {
        		if (includeBond[mConnBond[atom][i]]) {
        			atomMap[atom] = copyAtom(destMol, atom, 0, 0);
        			break;
        			}
        		}
        	}

        destMol.mAllBonds = 0;
        for (int bnd=0; bnd<mAllBonds;bnd++)
            if (includeBond[bnd]) {
                copyBond(destMol, bnd, 0, 0, atomMap, recognizeDelocalizedBonds);
            	}
	    	// if we keep one atom of a removed bi-polar bonds, then we need to neutralize
	    	else {
	        	int atom1 = mBondAtom[0][bnd];
	        	int atom2 = mBondAtom[1][bnd];
	        	if (atomMap[atom1] == -1
	        	  ^ atomMap[atom2] == -1) {
		    		if (mAtomCharge[atom1] != 0
		  			 && mAtomCharge[atom2] != 0
		    		 && ((mAtomCharge[atom1] < 0)
		    		   ^ (mAtomCharge[atom2] < 0))) {
			    		if (atomMap[atom1] != -1)
			    			destMol.mAtomCharge[atomMap[atom1]] += (mAtomCharge[atom1] < 0) ? 1 : -1;
			    		if (atomMap[atom2] != -1)
			    			destMol.mAtomCharge[atomMap[atom2]] += (mAtomCharge[atom2] < 0) ? 1 : -1;
		    			}
	        		}
	    		}

        copyMoleculeProperties(destMol);
        destMol.mValidHelperArrays = cHelperNone;

        destMol.renumberESRGroups(cESRTypeAnd);
        destMol.renumberESRGroups(cESRTypeOr);

        if (destMol.mAllAtoms != mAllAtoms)
            destMol.setFragment(true);

        if (recognizeDelocalizedBonds)
            new AromaticityResolver(destMol).locateDelocalizedDoubleBonds();
 
        return atomMap;
        }


	/**
	 * @param atom
	 * @return the number of direct neighbor atoms including explicit hydrogen atoms
	 */
	public int getAllConnAtoms(int atom) {
		return mAllConnAtoms[atom];
		}


	/**
	 * @param atom
	 * @return the number of connected explicit and implicit hydrogen atoms
	 */
	public int getAllHydrogens(int atom) {
		return mAllConnAtoms[atom] - mConnAtoms[atom] + getImplicitHydrogens(atom);
		}


	/**
	 * @return the number of non-hydrogen atoms in the molecule
	 */
	public int getAtoms() {
		return mAtoms;
		}


	/**
	 * This is different from the Hendrickson pi-value, which considers pi-bonds to carbons only.
	 * @param atom
	 * @return the number pi electrons at atom (the central atom of acetone would have 1)
	 */
	public int getAtomPi(int atom) {
		return mPi[atom];
		}


	/**
	 * @param atom
	 * @return Hendrickson sigma-value, which is the number attached carbon atoms
	 *
	public int getAtomSigma(int atom) {
		int sigma = 0;
		for (int i=0; i<mConnAtoms[atom]; i++)
			if (mAtomicNo[mConnAtom[atom][i]] == 6)
				sigma++;
		return sigma;
		}*/


	/**
	 * @param atom
	 * @return Hendrickson Z-value, which is the sum of all bond orders to any attached hetero atoms
	 *
	public int getAtomZValue(int atom) {
		int z = 0;
		for (int i=0; i<mConnAtoms[atom]; i++)
			if (isElectronegative(mConnAtom[atom][i]))
				z += mConnBondOrder[atom][i];
		return z;
		}*/

	
	/**
	 * @param atom
	 * @return 0 or the size of the smallest ring that atom is a member of
	 */
	public int getAtomRingSize(int atom) {
		return (mRingSet != null && atom<mAtoms) ?
                mRingSet.getAtomRingSize(atom) : 0;
		}


	/**
	 * @param bond
	 * @return 0 or the size of the smallest ring that bond is a member of
	 */
    public int getBondRingSize(int bond) {
        return (mRingSet != null && bond<mBonds) ?
                mRingSet.getBondRingSize(bond) : 0;
        }


	/**
	 * @return the number of bonds in the molecule not counting those that connect natural-abundance hydrogen atoms
	 */
	public int getBonds() {
		return mBonds;
		}


	/**
	 * @return -1 or the bond that connects atom1 with atom2
	 */
	public int getBond(int atom1, int atom2) {
		for (int i=0; i<mAllConnAtoms[atom1]; i++)
			if (mConnAtom[atom1][i] == atom2)
				return mConnBond[atom1][i];
		return -1;
		}


	/**
	 * @return a copy of this with all arrays sized to just cover all existing atoms and bonds
	 */
	public ExtendedMolecule getCompactCopy() {
		ExtendedMolecule theCopy = new ExtendedMolecule(mAllAtoms, mAllBonds);
		copyMolecule(theCopy);
		return theCopy;
		}


	/**
	 * @param atom
	 * @param i
	 * @return the i-th neighbor atom of atom
	 */
	public int getConnAtom(int atom, int i) {
		return mConnAtom[atom][i];
		}


	/**
	 * @param atom
	 * @return the number of direct neighbor atoms excluding hydrogen atoms
	 */
	public int getConnAtoms(int atom) {
		return mConnAtoms[atom];
		}


	/**
	 * @param atom
	 * @param i
	 * @return the bond connecting atom and its i-th neighbor atom
	 */
	public int getConnBond(int atom, int i) {
		return mConnBond[atom][i];
		}


	/**
	 * @param atom
	 * @param i
	 * @return the order of the bond connecting atom and its i-th neighbor atom
	 */
	public int getConnBondOrder(int atom, int i) {
		return mConnBondOrder[atom][i];
		}


	public float getAverageBondLength(boolean nonHydrogenBondsOnly) {
        if (nonHydrogenBondsOnly) {
        	ensureHelperArrays(cHelperNeighbours);
    		return getAverageBondLength(mAtoms, mBonds);
        	}
        else {
    		return getAverageBondLength(mAllAtoms, mAllBonds);
        	}
		}


	/**
	 * Creates an array that maps connAtoms/connBonds sorted by atom indices.
	 * getConnAtom(atom, getSortedConnMap(atom, 0)) retrieves that neighbour
	 * of atom with the lowest atom index, i.e. that is the first in the atom table.
	 * @return neighbour index map
	 */
	public int[] getSortedConnMap(int atom) {
        int[] indexMap = new int[mAllConnAtoms[atom]];
        for (int i=0; i<mAllConnAtoms[atom]; i++)
            indexMap[i] = (mConnAtom[atom][i] << 16) + i;
        java.util.Arrays.sort(indexMap);
        for (int i=0; i<mAllConnAtoms[atom]; i++)
            indexMap[i] &= 0x0000FFFF;
        return indexMap;
	    }


	/**
	 * The sum of bond orders of explicitly connected neighbour atoms including hydrogen.
     * @param atom
     * @return explicitly used valence
	 */
	public int getOccupiedValence(int atom) {
		ensureHelperArrays(cHelperNeighbours);

		int valence = 0;
		for (int i=0; i<mAllConnAtoms[atom]; i++)
			valence += mConnBondOrder[atom][i];

		return valence;
		}


	/**
	 * The free valence is the number of potential additional single bonded
	 * neighbours to reach the atom's maximum valence. Atom charges are considered.
	 * Implicit hydrogens are not considered.
	 * Thus, the oxygen in a R-O(-) has a free valence of 0, the nitrogen in R3N(+)
	 * has a free valence of 1.
	 * @param atom
	 * @return
	 */
	public int getFreeValence(int atom) {
		return getMaxValence(atom) - getOccupiedValence(atom);
		}


	/**
	 * If the explicitly attached neighbors cause an atom valence to exceed
	 * the lowest allowed valence for this atomic no, then this method returns
	 * the next higher allowed valence, e.g. O=P(-H)-OMe :<br>
	 * standard P valence is 3, used valence is 4, implicit abnormal valence is 5.
	 * The molecule is interpreted as O=PH2-OMe. Requires cHelperNeighbours!
	 * @param atom
	 * @param neglectExplicitHydrogen
	 * @return abnormal valence or -1 if valence doesn't exceed standard valence
	 */
	public int getImplicitHigherValence(int atom, boolean neglectExplicitHydrogen) {
	    int occupiedValence = getOccupiedValence(atom) - getElectronValenceCorrection(atom);
	    if (neglectExplicitHydrogen)
	        occupiedValence -= mAllConnAtoms[atom] - mConnAtoms[atom];

	    byte[] valenceList = (mAtomicNo[atom] < cAtomValence.length) ?
                cAtomValence[mAtomicNo[atom]] : null;

        int valence = (valenceList == null) ? cDefaultAtomValence : valenceList[0];
        if (occupiedValence <= valence)
            return -1;

        // If we don't neglect hydrogen and don't have allowed higher valences
        // then consider explicit hydrogens as errors.
        if (!neglectExplicitHydrogen
         && (valenceList == null || valenceList.length == 1)) {
            occupiedValence -= mAllConnAtoms[atom] - mConnAtoms[atom];
            return (occupiedValence <= valence) ? -1 : occupiedValence;
            }

        if (valenceList != null)
        	for (int i=1; (valence<occupiedValence) && (i<valenceList.length); i++)
        		valence = valenceList[i];

        // If we don't neglect hydrogen and the maximum allowed higher valence
        // is still exceeded then consider explicit hydrogens as errors.
        if (!neglectExplicitHydrogen)
            occupiedValence -= mAllConnAtoms[atom] - mConnAtoms[atom];

        return Math.max(valence, occupiedValence);
	    }


	/**
	 * Calculates for every atom the mean value of all shortest routes (bonds in between)
	 * to any other atom of the same fragment.
	 * @return 
	 */
    public float[] getAverageTopologicalAtomDistance() {
        ensureHelperArrays(cHelperNeighbours);

        float[] meanDistance = new float[mAtoms];
		int[] graphAtom = new int[mAtoms];
		for (int startAtom=0; startAtom<mAtoms; startAtom++) {
			graphAtom[0] = startAtom;
			int[] graphLevel = new int[mAtoms];
			graphLevel[startAtom] = 1;

			int current = 0;
			int highest = 0;
			while (current <= highest) {
				for (int i=0; i<mConnAtoms[graphAtom[current]]; i++) {
					int candidate = mConnAtom[graphAtom[current]][i];
					if (graphLevel[candidate] == 0) {
						graphLevel[candidate] = graphLevel[graphAtom[current]] + 1;
						graphAtom[++highest] = candidate;
						meanDistance[startAtom] += (graphLevel[candidate] - 1);
						}
					}
				current++;
				}
			meanDistance[startAtom] /= highest;
			}

		return meanDistance;
    	}

    /**
     * Calculates the length of the shortest path between atoms atom1 and atom2
     * @param atom1
     * @param atom2
     * @return path length (no of bonds); -1 if there is no path
     */
    public int getPathLength(int atom1, int atom2) {
        if (atom1 == atom2)
            return 0;

        ensureHelperArrays(cHelperNeighbours);

        int[] graphLevel = new int[mAtoms];
        int graphAtom[] = new int[mAtoms];

        graphAtom[0] = atom1;
        graphLevel[atom1] = 1;
        int current = 0;
        int highest = 0;
        while (current <= highest) {
            for (int i=0; i<mConnAtoms[graphAtom[current]]; i++) {
                int candidate = mConnAtom[graphAtom[current]][i];
                if (candidate == atom2)
                    return graphLevel[graphAtom[current]];
                if (graphLevel[candidate] == 0) {
                    graphAtom[++highest] = candidate;
                    graphLevel[candidate] = graphLevel[graphAtom[current]]+1;
                    }
                }
            current++;
            }
        return -1;
        }


    /**
     * Calculates the length of the shortest path between atoms atom1 and atom2,
     * which is not larger than maxLength and avoids atoms indicated by neglectAtom.
     * @param atom1
     * @param atom2
     * @param maxLength paths larger than maxLength won't be detected
     * @param neglectAtom null or atoms flagged which are forbidden path members
     * @return path length (no of bonds); -1 if there is no path
     */
    public int getPathLength(int atom1, int atom2, int maxLength, boolean[] neglectAtom) {
        if (atom1 == atom2)
            return 0;

        ensureHelperArrays(cHelperNeighbours);

        int[] graphLevel = new int[mAtoms];
        int graphAtom[] = new int[mAtoms];

        graphAtom[0] = atom1;
        graphLevel[atom1] = 1;
        int current = 0;
        int highest = 0;
        while (current <= highest && graphLevel[graphAtom[current]] <= maxLength) {
            for (int i=0; i<mConnAtoms[graphAtom[current]]; i++) {
                int candidate = mConnAtom[graphAtom[current]][i];
                if (candidate == atom2)
                    return graphLevel[graphAtom[current]];

                if (graphLevel[candidate] == 0 && (neglectAtom == null || !neglectAtom[candidate])) {
                    graphAtom[++highest] = candidate;
                    graphLevel[candidate] = graphLevel[graphAtom[current]]+1;
                    }
                }
            current++;
            }
        return -1;
        }


    /**
     * Locates and returns the shortest path between atoms atom1 and atom2
     * @param pathAtom array large enough to hold all path atoms, i.e. maxLength+1
     * @param atom1 first atom of path; ends up in pathAtom[0]
     * @param atom2 last atom of path; ends up in pathAtom[pathLength]
     * @param maxLength paths larger than maxLength won't be detected
     * @param neglectBond null or bitmask of forbidden bonds
     * @return number of bonds of path; -1 if there is no path
     */
    public int getPath(int[] pathAtom, int atom1, int atom2, int maxLength, boolean[] neglectBond) {
        if (atom1 == atom2) {
            pathAtom[0] = atom1;
            return 0;
            }

        ensureHelperArrays(cHelperNeighbours);

        int[] graphLevel = new int[mAtoms];
        int graphAtom[] = new int[mAtoms];
        int parentAtom[] = new int[mAtoms];

        graphAtom[0] = atom1;
        graphLevel[atom1] = 1;

        int current = 0;
        int highest = 0;
        while (current <= highest && graphLevel[graphAtom[current]] <= maxLength) {
            int parent = graphAtom[current];
            for (int i=0; i<mConnAtoms[parent]; i++) {
                if (neglectBond == null || !neglectBond[mConnBond[parent][i]]) {
                    int candidate = mConnAtom[parent][i];
                    if (candidate == atom2) {
                        int index = graphLevel[parent];
                        pathAtom[index] = candidate;
                        pathAtom[--index] = parent;
                        while (index > 0) {
                            pathAtom[index-1] = parentAtom[pathAtom[index]];
                            index--;
                            }
                        return graphLevel[parent];
                        }
    
                    if (graphLevel[candidate] == 0) {
                        graphAtom[++highest] = candidate;
                        graphLevel[candidate] = graphLevel[parent]+1;
                        parentAtom[candidate] = parent;
                        }
                    }
                }
            current++;
            }
        return -1;
        }


    /**
     * Finds bonds of a path that is defined by an atom sequence.
     * @param pathAtom pathAtom[0]...[pathLength] -> list of atoms on path 
     * @param pathBond int array not smaller than pathLength
     * @param pathLength no of path bonds == no of path atoms - 1
     */
    public void getPathBonds(int[] pathAtom, int[] pathBond, int pathLength) {
        ensureHelperArrays(cHelperNeighbours);
        for (int i=0; i<pathLength; i++) {
            for (int j=0; j<mConnAtoms[pathAtom[i]]; j++) {
                if (mConnAtom[pathAtom[i]][j] == pathAtom[i+1]) {
                    pathBond[i] = mConnBond[pathAtom[i]][j];
                    break;
                    }
                }
            }
        }


    /**
     * @param atom1
     * @param atom2
     * @return whether there is a path of bonds leading from atom1 to atom2
     */
    public boolean shareSameFragment(int atom1, int atom2) {
		return (getPathLength(atom1, atom2) != -1);
		}


    /**
     * This adds a fragment from sourceMol to this molecule by first copying rootAtom and then
     * all connected atoms and bonds by traversing the graph breadth first.
     * @param sourceMol molecule from which the fragment is copied to this
     * @param rootAtom
	 * @param atomMap null or int[] not smaller than sourceMol.mAllAtoms; receives atom indices of this molecule
     */
	public void addFragment(ExtendedMolecule sourceMol, int rootAtom, int[] atomMap) {
		sourceMol.ensureHelperArrays(cHelperNeighbours);

        if (atomMap == null)
        	atomMap = new int[sourceMol.mAllAtoms];

        int esrGroupCountAND = renumberESRGroups(Molecule.cESRTypeAnd);
        int esrGroupCountOR = renumberESRGroups(Molecule.cESRTypeOr);

		boolean[] isFragmentMember = new boolean[sourceMol.mAllAtoms];
		int graphAtom[] = new int[sourceMol.mAllAtoms];

		graphAtom[0] = rootAtom;
		isFragmentMember[rootAtom] = true;
		atomMap[rootAtom] = sourceMol.copyAtom(this, rootAtom, esrGroupCountAND, esrGroupCountOR);

		int current = 0;
		int highest = 0;
	 	while (current <= highest) {
			for (int i=0; i<sourceMol.mAllConnAtoms[graphAtom[current]]; i++) {
				int candidate = sourceMol.mConnAtom[graphAtom[current]][i];
				if (!isFragmentMember[candidate]) {
					graphAtom[++highest] = candidate;
					isFragmentMember[candidate] = true;
					atomMap[candidate] = sourceMol.copyAtom(this, candidate, esrGroupCountAND, esrGroupCountOR);
					}
				}
			current++;
			}

        for (int bond=0; bond<sourceMol.mAllBonds; bond++)
           	if (isFragmentMember[sourceMol.mBondAtom[0][bond]])
           		sourceMol.copyBond(this, bond, esrGroupCountAND, esrGroupCountOR, atomMap, false);

        renumberESRGroups(cESRTypeAnd);
        renumberESRGroups(cESRTypeOr);

        mValidHelperArrays = cHelperNone;
		}


	/**
     * @param rootAtom
     * @return array of all atoms that for which a path of bonds leads to rootAtom
     */
	public int[] getFragmentAtoms(int rootAtom) {
		ensureHelperArrays(cHelperNeighbours);

		boolean[] isFragmentMember = new boolean[mAllAtoms];
		int graphAtom[] = new int[mAllAtoms];

		graphAtom[0] = rootAtom;
		isFragmentMember[rootAtom] = true;
		int current = 0;
		int highest = 0;
		int fragmentMembers = 1;
	 	while (current <= highest) {
			for (int i=0; i<mAllConnAtoms[graphAtom[current]]; i++) {
				int candidate = mConnAtom[graphAtom[current]][i];
				if (!isFragmentMember[candidate]) {
					graphAtom[++highest] = candidate;
					isFragmentMember[candidate] = true;
					fragmentMembers++;
					}
				}
			current++;
			}

		int[] fragmentMember = new int[fragmentMembers];
		fragmentMembers = 0;
		for (int atom=0; atom<mAllAtoms; atom++)
			if (isFragmentMember[atom])
				fragmentMember[fragmentMembers++] = atom;

		return fragmentMember;
		}


	/**
	 * Locates all disconnected fragments in the molecule and assigns
	 * fragment numbers (starting from 0) to all atoms. Individual bonds can be marked to
	 * be skipped, i.e. to be treated as non-existing.
	 * @param fragmentNo array not smaller than getAllAtoms()
	 * @param neglectBond array not smaller than getAllBonds()
	 * @return array of atom's fragments indexes
	 */
    public int getFragmentNumbers(int[] fragmentNo, boolean[] neglectBond) {
        ensureHelperArrays(cHelperNeighbours);

        for (int atom=0; atom<mAllAtoms; atom++)
            fragmentNo[atom] = -1;

        int fragments = 0;
        for (int atom=0; atom<mAllAtoms; atom++) {
            if (fragmentNo[atom] == -1) {
                fragmentNo[atom] = fragments;
                int graphAtom[] = new int[mAllAtoms];
                graphAtom[0] = atom;
                int current = 0;
                int highest = 0;
                while (current <= highest) {
                    for (int i=0; i<mAllConnAtoms[graphAtom[current]]; i++) {
                        int candidate = mConnAtom[graphAtom[current]][i];
                        if (fragmentNo[candidate] == -1
                         && !neglectBond[mConnBond[graphAtom[current]][i]]) {
                            graphAtom[++highest] = candidate;
                            fragmentNo[candidate] = fragments;
                            }
                        }
                    current++;
                    }
                fragments++;
                }
            }
        return fragments;
        }


    /**
     * Locates all unconnected fragments in the Molecule and assigns fragment indexes
     * for every atom starting with 0. Optionally the fragment detection may be restricted to
     * those atoms that have been previously marked with setAtomMarker(). In that case
     * non-marked atoms receive the fragment number -1 and are not considered a connection between
     * marked atoms potentially causing two marked atoms to end up in different fragments, despite
     * sharing the same fragment.
     * @param fragmentNo array at least mAllAtoms big to receive atom fragment indexes
     * @param markedAtomsOnly if true, then only atoms marked with setAtomMarker() are considered
     * @return number of disconnected fragments
     */
    public int getFragmentNumbers(int[] fragmentNo, boolean markedAtomsOnly) {
		ensureHelperArrays(cHelperNeighbours);

		for (int atom=0; atom<mAllAtoms; atom++)
			fragmentNo[atom] = -1;

		int fragments = 0;
		for (int atom=0; atom<mAllAtoms; atom++) {
			if (fragmentNo[atom] == -1
             && (!markedAtomsOnly || isMarkedAtom(atom))) {
				fragmentNo[atom] = fragments;
				int graphAtom[] = new int[mAllAtoms];
				graphAtom[0] = atom;
				int current = 0;
				int highest = 0;
				while (current <= highest) {
					for (int i=0; i<mAllConnAtoms[graphAtom[current]]; i++) {
						int candidate = mConnAtom[graphAtom[current]][i];
						if (fragmentNo[candidate] == -1
                         && (!markedAtomsOnly || isMarkedAtom(candidate))) {
							graphAtom[++highest] = candidate;
							fragmentNo[candidate] = fragments;
							}
						}
					current++;
					}
				fragments++;
				}
			}
		return fragments;
		}


	/**
	 * Removes all unconnected fragments except for the largest one.
	 * If small fragments were removed, then canonizeCharge() is called to
	 * neutralize charges after potential removal of counter ions.
	 * @return true, if fragments were removed
	 */
    public boolean stripSmallFragments() {
    	int[] fragmentNo = new int[mAllAtoms];
    	int fragmentCount = getFragmentNumbers(fragmentNo, false);

		if (fragmentCount <= 1)
			return false;

		int[] fragmentSize = new int[fragmentCount];
		for (int atom=0; atom<mAllAtoms; atom++)
			fragmentSize[fragmentNo[atom]]++;

		int largestFragment = 0;
		int largestSize = fragmentSize[0];
		for (int i=1; i<fragmentCount; i++) {
			if (largestSize < fragmentSize[i]) {
				largestSize = fragmentSize[i];
				largestFragment = i;
				}
			}

		for (int atom=0; atom<mAllAtoms; atom++)
			if (fragmentNo[atom] != largestFragment)
				mAtomicNo[atom] = -1;				// mark for delete

		for (int bond=0; bond<mAllBonds; bond++)
			if (fragmentNo[mBondAtom[0][bond]] != largestFragment)
				mBondType[bond] = cBondTypeDeleted;	// mark for delete

		compressMolTable();
		mValidHelperArrays = cHelperNone;

		try { canonizeCharge(true); } catch (Exception e) {}

		return true;
		}


    /**
	 * Starting from startAtom this method locates a system of annelated or bridged ring systems
	 * with all members bonds being a ring bond. Detected member atoms and bonds are flagged
	 * accordingly.
	 * @param startAtom
	 * @param aromaticOnly if set then only aromatic atoms and bonds are considered
	 * @param isMemberAtom
	 * @param isMemberBond
	 */
	public void findRingSystem(int startAtom, boolean aromaticOnly, boolean[] isMemberAtom, boolean[] isMemberBond) {
		ensureHelperArrays(cHelperRings);
		
		if (!isRingAtom(startAtom) || (aromaticOnly && !isAromaticAtom(startAtom)))
			return;

		int[] graphAtom = new int[mAtoms];
		    
		graphAtom[0] = startAtom;
		isMemberAtom[startAtom] = true;

		int current = 0;
		int highest = 0;
	 	while (current <= highest) {
			for (int i=0; i<mConnAtoms[graphAtom[current]]; i++) {
				int candidateBond = mConnBond[graphAtom[current]][i];
				if (!isMemberBond[candidateBond]
				 && isRingBond(candidateBond)
				 && (!aromaticOnly || isAromaticBond(candidateBond))) {
					isMemberBond[candidateBond] = true;
					int candidateAtom = mConnAtom[graphAtom[current]][i];
					if (!isMemberAtom[candidateAtom]) {
					    isMemberAtom[candidateAtom] = true;
						graphAtom[++highest] = candidateAtom;
						}
					}
				}
			current++;
			}
		}

	/**
     * Determines all atoms of the substituent attached to coreAtom and starting
     * with firstAtom. If isMemberAtom!=null, then all substituent member atoms
     * will have the the respective index being flagged upon return. This includes
     * firstAtom and excludes coreAtom.
     * If substituent!=null, then it will contain the substituent as Molecule.
     * At the position of the coreAtom substituent will contain a wildcard atom.
	 * If substituent!=null and atomMap!=null then atomMap receives atom index mapping from
	 * this to substituent with non-member atoms being -1.
     * Returns -1 and an empty substituent if coreAtom and firstAtom share a ring
     * @param coreAtom the atom to which the substituent is connected
     * @param firstAtom the substituent's atom that is connected to coreAtom
     * @param isMemberAtom may be null, otherwise set to contain atom membership mask
     * @param substituent may be null, otherwise set to contain the substituent
     * @param atomMap null or int[] not smaller than this.getAllAtoms()
     * @return substituent atom count not counting coreAtom; -1 if coreAtom and firstAtom share a ring
     */
	public int getSubstituent(int coreAtom, int firstAtom, boolean[] isMemberAtom, ExtendedMolecule substituent, int[] atomMap) {
		ensureHelperArrays(cHelperNeighbours);

		if (substituent != null) {
    		substituent.deleteMolecule();
    		substituent.mIsFragment = false;
		    }

		int[] graphAtom = new int[mAllAtoms];
		if (isMemberAtom == null)
		    isMemberAtom = new boolean[mAllAtoms];
		else
		    java.util.Arrays.fill(isMemberAtom, false);
		    
		graphAtom[0] = coreAtom;
		graphAtom[1] = firstAtom;
		isMemberAtom[coreAtom] = true;
		isMemberAtom[firstAtom] = true;
		int current = 1;
		int highest = 1;
	 	while (current <= highest) {
			for (int i=0; i<mAllConnAtoms[graphAtom[current]]; i++) {
				int candidate = mConnAtom[graphAtom[current]][i];
				if (candidate == coreAtom) {
					if (current != 1)
						return -1;
					}
				if (!isMemberAtom[candidate]) {
				    isMemberAtom[candidate] = true;
					graphAtom[++highest] = candidate;
					}
				}
			current++;
			}

	 	if (substituent != null) {
	 		if (atomMap == null)
	 			atomMap = new int[isMemberAtom.length];
	 	    copyMoleculeByAtoms(substituent, isMemberAtom, false, atomMap);
	        substituent.changeAtom(atomMap[coreAtom], 0, 0, -1, 0);
	 	    }

        isMemberAtom[coreAtom] = false;
		return highest;
		}

	/**
	 * Counts the number of atoms of the substituent connected to coreAtom
	 * defined by firstAtom and not including the coreAtom.
	 * @param coreAtom
	 * @param firstAtom
	 * @return atom count of substituent or -1 if coreAtom and firstAtom are in the same ring
	 */
    public int getSubstituentSize(int coreAtom, int firstAtom) {
        ensureHelperArrays(cHelperNeighbours);
    
        int[] graphAtom = new int[mAtoms];
        boolean[] isMember = new boolean[mAtoms];
        graphAtom[0] = coreAtom;
        graphAtom[1] = firstAtom;
        isMember[coreAtom] = true;
        isMember[firstAtom] = true;
        int current = 1;
        int highest = 1;
        while (current <= highest) {
            for (int i=0; i<mConnAtoms[graphAtom[current]]; i++) {
                int candidate = mConnAtom[graphAtom[current]][i];
                if (candidate == coreAtom) {
                    if (current != 1)
                        return -1;
                    }
                if (!isMember[candidate]) {
                    isMember[candidate] = true;
                    graphAtom[++highest] = candidate;
                    }
                }
            current++;
            }
        return highest;
        }

    /**
     * 
     * @param atom
     * @return
     */
	public int getImplicitHydrogens(int atom) {
		if (mIsFragment
		 && (mAtomQueryFeatures[atom] & cAtomQFNoMoreNeighbours) == 0)
				return 0;

		int atomicNo = mAtomicNo[atom];
		if (atomicNo != 5    // B
		 && atomicNo != 6    // C
         && atomicNo != 7    // N
         && atomicNo != 8    // O
         && atomicNo != 9    // F
         && atomicNo != 13   // Al
         && atomicNo != 14   // Si
         && atomicNo != 15   // P
         && atomicNo != 16   // S
         && atomicNo != 17   // Cl
         && atomicNo != 33   // As
         && atomicNo != 34   // Se
         && atomicNo != 35   // Br
         && atomicNo != 52   // Te
         && atomicNo != 53)  // I
		    return 0;

		ensureHelperArrays(cHelperNeighbours);

        int occupiedValence = 0;
		for (int i=0; i<mAllConnAtoms[atom]; i++)
		    occupiedValence += mConnBondOrder[atom][i];

		if (mIsFragment) {
			int delocalizedBonds = 1;
			for (int i=0; i<mConnAtoms[atom]; i++)
				if (mBondType[mConnBond[atom][i]] == cBondTypeDelocalized)
					delocalizedBonds++;
			occupiedValence += delocalizedBonds >> 1;
			}

		occupiedValence -= getElectronValenceCorrection(atom);
        int maxValence = getAtomAbnormalValence(atom);
        if (maxValence == -1) {
            if (mAtomicNo[atom] >= 171 && mAtomicNo[atom] <= 190) {
                maxValence = 2;
                }
            else {
                byte[] valenceList = (mAtomicNo[atom] < cAtomValence.length) ?
                        cAtomValence[mAtomicNo[atom]] : null;
                if (valenceList == null) {
                    maxValence = cDefaultAtomValence;
                    }
                else {
                    maxValence = valenceList[0];
                    for (int i=1; (maxValence<occupiedValence) && (i<valenceList.length); i++)
                        maxValence = valenceList[i];
                    }
                }
            }

		return Math.max(0, maxValence - occupiedValence);
		}

	/**
	 * Calculates a rounded mass of the molecule
	 * @return
	 */
	public int getMolweight() {
		ensureHelperArrays(cHelperNeighbours);
		int molweight = 0;
		for (int atom=0; atom<mAllAtoms; atom++) {
			int mass = mAtomMass[atom] != 0 ? mAtomMass[atom] : cRoundedMass[mAtomicNo[atom]];
            molweight += mass + getImplicitHydrogens(atom) * cRoundedMass[1];
            if (mAtomicNo[atom] >= 171 && mAtomicNo[atom] <= 190
			 && mAllConnAtoms[atom] > 2)	// to accomodate for bisulfide bridges
				molweight -= (mAllConnAtoms[atom] - 2) * cRoundedMass[1];
			}

        return molweight;
        }

	/**
	 * Simple method to calculate rotatable bonds. This method counts all single
	 * bonds provided that they<br>
	 * - are not a terminal bond<br>
	 * - are not part of a ring<br>
	 * - are not an amide bond<br>
	 * - are not the second of two equivalent bonds next to the same triple bond<br>
	 * @return
	 */
    public int getRotatableBondCount() {
        int rCount = 0;
        ensureHelperArrays(Molecule.cHelperRings);
        for (int bond=0; bond<mBonds; bond++) {
            if (getBondOrder(bond) == 3
             && mConnAtoms[mBondAtom[0][bond]] > 1
             && mConnAtoms[mBondAtom[1][bond]] > 1) {
                rCount++;   // count x-x#x-x as one rotatable bond
                continue;
                }

            if (getBondOrder(bond) != 1 || isRingBond(bond))
                continue;

            boolean isRotatable = true;
            for (int i=0; i<2; i++) {
                int atom1 = mBondAtom[i][bond];
                if (mConnAtoms[atom1] == 1) {
                    isRotatable = false;
                    break;  // terminal bond
                    }
                if (mPi[atom1] == 2
                 && mConnAtoms[atom1] == 2) {
                    isRotatable = false;
                    break;  // vicinal to triple bond (count triple bond instead)
                    }
                if (mAtomicNo[atom1] == 7 && !isAromaticAtom(atom1)) {
                    int atom2 = mBondAtom[1-i][bond];
                    for (int j=0; j<mConnAtoms[atom2]; j++) {
                        int connAtom = mConnAtom[atom2][j];
                        int connBond = mConnBond[atom2][j];
                        if (connBond != bond
                         && getBondOrder(connBond) > 1
                         && !isAromaticAtom(connAtom)
                         && isElectronegative(connAtom)) {
                            isRotatable = false;
                            break;  // amid bond
                            }
                        }
                    }
                }

            if (isRotatable)
                rCount++;
            }
        return rCount;
        }

    public int getAromaticRingCount() {
        ensureHelperArrays(cHelperRings);
    	int count = 0;
    	for (int i=0; i<mRingSet.getSize(); i++)
    		if (mRingSet.isAromatic(i))
    			count++;
    	return count;
    	}

    /**
     * Calculates the number of independent rings of which 'atom' is a member.
     * Any combination of two connected atoms to 'atom' is used for:
     * - finding the shortest path connecting these two neighbors avoiding 'atom'
     * - if such a path exists and at least one bonds of that path is not a member
     *   of a path found earlier then count this path as an independent ring closure.
     * @param atom
     * @param maxRingSize
     * @return number of independent rings
     */
    public int getAtomRingCount(int atom, int maxRingSize) {
        ensureHelperArrays(cHelperRings);
        boolean[] bondTouched = new boolean[mBonds];
        boolean[] neglectBond = new boolean[mBonds];
        int[] ringAtom = new int[mAtoms];
        int count = 0;
        for (int i=1; i<mConnAtoms[atom]; i++) {
            int bond1 = mConnBond[atom][i];
            if (isRingBond(bond1)) {
                for (int j=0; j<i; j++) {
                    int bond2 = mConnBond[atom][j];
                    if (isRingBond(bond2)) {
                        neglectBond[bond1] = true;
                        neglectBond[bond2] = true;
                        int pathLength = getPath(ringAtom, mConnAtom[atom][i], mConnAtom[atom][j], maxRingSize-2, neglectBond);
                        neglectBond[bond1] = false;
                        neglectBond[bond2] = false;
                        if (pathLength != -1) {
                            boolean isIndependentRing = false;
                            int[] pathBond = new int[pathLength];
                            getPathBonds(ringAtom, pathBond, pathLength);
                            for (int k=0; k<pathLength; k++) {
                                if (!bondTouched[pathBond[k]]) {
                                    bondTouched[pathBond[k]] = true;
                                    isIndependentRing = true;
                                    }
                                }
                            if (isIndependentRing)
                                count++;
                            }
                        }
                    }
                }
            }
        return count;
        }

    /**
     * @return a RingCollection object, which contains a total set of small rings
     */
    public RingCollection getRingSet() {
		ensureHelperArrays(cHelperRings);
		return mRingSet;
		}

    /**
     * Locates that single bond which is the preferred one to be converted into up/down bond
     * in order to define the atom chirality.
     * @param atom parity carrying atom, i.e. a tetrahedral stereocenter or central allene atom
     * @return preferred bond or -1, if no single bond existing
     */
	public int getAtomPreferredStereoBond(int atom) {
		ensureHelperArrays(cHelperRings);
        if (mPi[atom] == 2 && mConnAtoms[atom] == 2)
            return preferredAlleneStereoBond(atom);
        else
			return preferredTHStereoBond(atom);
		}


    /**
     * Locates that single bond which is the preferred one to be converted into up/down bond
     * in order to define the bond chirality.
     * @param bond BINAP type of chirality bond
     * @return preferred bond or -1, if no single bond existing
     */
	public int getBondPreferredStereoBond(int bond) {
		return preferredBinapStereoBond(bond);
		}


	private int getStereoBondScore(int bond, int atom) {
			// score used to select one bond to be the stereo bond
        if (getBondOrder(bond) != 1)
            return 0;

        return 16 - mAllConnAtoms[atom]
                  + ((mAtomicNo[atom] == 1) ? 4096 : 0)
                  + (((mBondType[bond] & cBondTypeMaskStereo) == 0 || mBondAtom[0][bond] != atom) ? 2048 : 0)
                  + ((getAtomParity(atom) == 0) ? 1024 : 0)
                  + ((!isRingBond(bond)) ? 512 : 0)
                  + ((mAtomicNo[atom] != 6) ? 256 : 0);
		}


	/**
	 * @param atom
	 * @return whether the atom is in an allylic/benzylic position
	 */
	public boolean isAllylicAtom(int atom) {
		return (mAtomFlags[atom] & cAtomFlagAllylic) != 0;
		}


	public boolean isAromaticAtom(int atom) {
		return (mAtomFlags[atom] & cAtomFlagAromatic) != 0;
		}


	public boolean isAromaticBond(int bnd) {
		return (mBondFlags[bnd] & cBondFlagAromatic) != 0;
		}


	/**
	 * A bond is considered delocalized, if it has different bond orders in
	 * different, but energetically equivalent mesomeric structures. Bonds in aromatic 6-membered
	 * rings typically are delocalized, while those in uncharged 5-membered aromatic rings are not.
	 * Indole has 6 delocalized bonds.
	 * @param bond
	 * @return
	 */
	public boolean isDelocalizedBond(int bond) {
		return (mBondFlags[bond] & cBondFlagDelocalized) != 0;
		}


	public boolean isRingAtom(int atom) {
		return (mAtomFlags[atom] & cAtomFlagsRingBonds) != 0;
		}


	public boolean isRingBond(int bnd) {
		return (mBondFlags[bnd] & cBondFlagRing) != 0;
		}


	/**
	 * @param atom
	 * @return whether atom is a member of a ring not larger than 7 atoms
	 */
	public boolean isSmallRingAtom(int atom) {
		return (mAtomFlags[atom] & cAtomFlagSmallRing) != 0;
		}


	/**
	 * @param bond
	 * @return whether bond is a member of a ring not larger than 7 atoms
	 */
	public boolean isSmallRingBond(int bond) {
		return (mBondFlags[bond] & cBondFlagSmallRing) != 0;
		}


	/**
	 * @param atom
	 * @return whether atom has a neighbor that is connected through a double/triple bond to a hetero atom
	 */
	public boolean isStabilizedAtom(int atom) {
		return (mAtomFlags[atom] & cAtomFlagStabilized) != 0;
		}


	public int getAtomRingBondCount(int atom) {
		int flags = (mAtomFlags[atom] & cAtomFlagsRingBonds);
		return (flags == 0) ? 0
			 : (flags == cAtomFlags2RingBonds) ? 2
			 : (flags == cAtomFlags3RingBonds) ? 3 : 4;
		}


	public String getChiralText() {
		return null;
		}


	/**
	 * Checks whether at least one of the connected bonds is a stereo bond.
	 * If atom is the central atom of an allene, then its direct neighbours
	 * are checked, whether one of them has a stereo bond.
	 * @param atom
	 * @return the stereo bond or -1 if not found
	 */
    public int getStereoBond(int atom) {
        ensureHelperArrays(cHelperNeighbours);
        if (mConnAtoms[atom] == 2
         && mConnBondOrder[atom][0] == 2
         && mConnBondOrder[atom][1] == 2) {
            for (int i=0; i<2; i++)
                for (int j=0; j<mAllConnAtoms[mConnAtom[atom][i]]; j++)
                    if (isStereoBond(mConnBond[mConnAtom[atom][i]][j], mConnAtom[atom][i]))
                        return mConnBond[mConnAtom[atom][i]][j];
            }
        else {
            for (int i=0; i<mAllConnAtoms[atom]; i++)
                if (isStereoBond(mConnBond[atom][i], atom))
                	return mConnBond[atom][i];
            }

        return -1;
        }


    /**
     * This tells the molecule that earlier calculated atom/bond parities are still valid.
     * If helperStereoBits!=0 then CIP parities and/or symmetry ranks can be declared valid as well.
     * Should be called if no coordinates are available but the parities are valid nevertheless.
     * Is used by the IDCodeParser after idcode parsing without coordinates.
     * (Note: After idcode parsing unknown stereo centers have parities cAtomParityNone
     * instead of cAtomParityUnknown and caling isStereoCenter(atom) returns false!!!)
     * This prevents the Canonizer to run the stereo recognition again when
     * ensureHelperArrays(cHelperParities or higher) is called.
     * May also be called after filling valences with hydrogen atoms, which have no coordinates,
     * to tell the molecule that the earlier created stereo flags are still valid.
	 * @param helperStereoBits 0 or combinations of cHelperBitCIP,cHelperBitSymmetry...,cHelperBitIncludeNitrogenParities
     */
    public void setParitiesValid(int helperStereoBits) {
        mValidHelperArrays |= (cHelperBitsStereo & (cHelperBitParities | helperStereoBits));
        }


    /**
     * This converts one single bond per parity into a stereo up/down bond to
     * correctly reflect the given parity. This works for tetrahedral and
     * allene atom parities as well as for BINAP type of bond parities.
     * Should only be called with valid TH and EZ parities and valid coordinates,
     * e.g. after idcode parsing with coordinates or after coordinate generation.
     */
	public void setStereoBondsFromParity() {
		ensureHelperArrays(cHelperRings); // in case we miss ring and neighbour information
    
		for (int atom=0; atom<mAtoms; atom++)
			setStereoBondFromAtomParity(atom);

		for (int bond=0; bond<mBonds; bond++)
			setStereoBondFromBondParity(bond);
		}


	/**
	 * Converts any stereo bond attached with its pointed tip
	 * to this atom into a single bond.
	 * @param atom
	 */
	public void convertStereoBondsToSingleBonds(int atom) {
        for (int i=0; i<mAllConnAtoms[atom]; i++) {
            int connBond = mConnBond[atom][i];
            if (isStereoBond(connBond, atom))
                mBondType[connBond] = cBondTypeSingle;
            }
	    }


	public void setStereoBondFromAtomParity(int atom) {
			// set an optimal bond to up/down to reflect the atom parity
		if (getAtomParity(atom) == Molecule.cAtomParityNone
		 || getAtomParity(atom) == Molecule.cAtomParityUnknown)
			return;

		if (mPi[atom] == 2 && mConnAtoms[atom] == 2) {
			setAlleneStereoBondFromParity(atom);
			return;
			}

        if (mConnAtoms[atom] < 3 || mConnAtoms[atom] > 4) {
            setAtomParity(atom, cAtomParityNone, false);
            return;
            }

        int[] sortedConnMap = getSortedConnMap(atom);

        float angle[] = new float[mAllConnAtoms[atom]];
        for (int i=0; i<mAllConnAtoms[atom]; i++)
            angle[i] = getBondAngle(mConnAtom[atom][sortedConnMap[i]], atom);

        for (int i=0; i<mAllConnAtoms[atom]; i++)
            if (mBondAtom[0][mConnBond[atom][i]] == atom
             && getBondOrder(mConnBond[atom][i]) == 1)
                mBondType[mConnBond[atom][i]] = cBondTypeSingle;

        if (setFisherProjectionStereoBondsFromParity(atom, sortedConnMap, angle))
            return;

        // If there is exactly one stereo bond at the atom then take this
        // as the new stereo bond. Otherwise find a preferred stereo bond.
        int preferredBond = -1;
        for (int i=0; i<mAllConnAtoms[atom]; i++) {
            int connBond = mConnBond[atom][i];
            if (isStereoBond(connBond, atom)) {
                mBondType[mConnBond[atom][i]] = cBondTypeSingle;
                if (preferredBond == -1)
                    preferredBond = connBond;
                else
                    preferredBond = -2;
                }
            }
        if (preferredBond < 0)
            preferredBond = preferredTHStereoBond(atom);

		if (mBondAtom[0][preferredBond] != atom) {
			mBondAtom[1][preferredBond] = mBondAtom[0][preferredBond];
			mBondAtom[0][preferredBond] = atom;
			}

		int preferredBondIndex = -1;
        for (int i=0; i<mAllConnAtoms[atom]; i++) {
            if (preferredBond == mConnBond[atom][sortedConnMap[i]]) {
                preferredBondIndex = i;
                break;
                }
            }

		final int[][] up_down = { { 2,1,2,1 },	// stereobond type to
		                          { 1,2,2,1 },	// achieve parity = 1
		                          { 1,1,2,2 },	// first dimension:
		                          { 2,1,1,2 },	// one of the 6 angle orders
		                          { 2,2,1,1 },	// second dimension:
		                          { 1,2,1,2 } };// index of ConnAtom bearing stereobond

/*		int preferredBondIndex = -1;
		int sortedConn[] = new int[4];
		for (int i=0; i<mAllConnAtoms[atom]; i++) {
			int lowestConnAtom = Integer.MAX_VALUE;
			int lowestConnIndex = -1;
			for (int j=0; j<mAllConnAtoms[atom]; j++) {
				if ((i == 0 || mConnAtom[atom][j] > sortedConn[i-1]) && lowestConnAtom > mConnAtom[atom][j]) {
					lowestConnAtom = mConnAtom[atom][j];
					lowestConnIndex = j;
					}
				}

			sortedConn[i] = lowestConnAtom;
			if (mConnBond[atom][lowestConnIndex] == preferredBond)
				preferredBondIndex = i;
			}    */

		for (int i=1; i<mAllConnAtoms[atom]; i++)
			if (angle[i] < angle[0])
				angle[i] += Math.PI*2;

		int bondType;
		if (mAllConnAtoms[atom] == 3) {
				// new handling!!! TLS 17.Oct.2005
				// Originally the parity depended solely on clockwise/anti-clockwise orientation
				// of three (priority sorted) angles, which included the angle of the stereo bond.
				// Now to handle strained rings with 'invalid' projections in an expected way
				// (all three bonds are within 180 degrees and the middle bond is the stereo bond
				// then treat this as if the stereo bond would be drawn 180 degrees rotated)
				// the angle of the stereobond is not considered anymore and the parity depends
				// solely on the question if the angle difference between higher priority bond
				// to lower priority bond is less or larger than 180 degrees.
			boolean inverted = false;
			switch (preferredBondIndex) {
			case 0:
				inverted = (((angle[1] < angle[2]) && (angle[2] - angle[1] < Math.PI))
						 || ((angle[1] > angle[2]) && (angle[1] - angle[2] > Math.PI)));
				break;
			case 1:
				inverted = (angle[2] - angle[0] > Math.PI);
				break;
			case 2:
				inverted = (angle[1] - angle[0] < Math.PI);
				break;
				}

			bondType = ((getAtomParity(atom) == cAtomParity1) ^ inverted) ?
						cBondTypeUp : cBondTypeDown;

/*				// original handling where parity depended solely on the clockwise/anti-clockwise of three angles
			bondType = ((getAtomParity(atom) == cAtomParity1)
					  ^ (angle[1] > angle[2])) ?
						cBondTypeDown : cBondTypeUp;	*/
			}
		else {
			int order = 0;
			if		(angle[1] <= angle[2] && angle[2] <= angle[3]) order = 0;
			else if (angle[1] <= angle[3] && angle[3] <= angle[2]) order = 1;
			else if (angle[2] <= angle[1] && angle[1] <= angle[3]) order = 2;
			else if (angle[2] <= angle[3] && angle[3] <= angle[1]) order = 3;
			else if (angle[3] <= angle[1] && angle[1] <= angle[2]) order = 4;
			else if (angle[3] <= angle[2] && angle[2] <= angle[1]) order = 5;
			bondType = ((getAtomParity(atom) == cAtomParity1)
					  ^ (up_down[order][preferredBondIndex] == 1)) ?
						cBondTypeDown : cBondTypeUp;
			}

		mBondType[preferredBond] = bondType;
		}


    private boolean setFisherProjectionStereoBondsFromParity(int atom, int[] sortedConnMap, float[] angle) {
        int[] direction = new int[mAllConnAtoms[atom]];
        int parity = getFisherProjectionParity(atom, sortedConnMap, angle, direction);
        if (parity == cAtomParityUnknown)
            return false;

        int bondType = (getAtomParity(atom) == parity) ? cBondTypeUp : cBondTypeDown;
        for (int i=0; i<mAllConnAtoms[atom]; i++) {
            if ((direction[i] & 1) == 1) {
                int bond = mConnBond[atom][sortedConnMap[i]];
                mBondType[bond] = bondType;
                if (mBondAtom[0][bond] != atom) {
                    mBondAtom[1][bond] = mBondAtom[0][bond];
                    mBondAtom[0][bond] = atom;
                    }
                }
            }

        return true;
        }


    /**
     * If the atom is a stereo center in fisher projection, then its
     * tetrahedral parity is returned. If the horizontal bonds are plain
     * single bonds, then they are interpreted as up-bonds.
     * @param atom the stereo center
     * @param index map of neighbours sorted by atom index
     * @param angle bond angles sorted by neighbour atom index
     * @param direction null or int[] large enough to receive bond directions
     * @return cAtomParity1,cAtomParity2 or cAtomParityUnknown
     */
    public int getFisherProjectionParity(int atom, int[] sortedConnMap, float[] angle, int[] direction) {
        if (direction == null)
            direction = new int[mAllConnAtoms[atom]];
        if (!getFisherProjectionBondDirections(atom, sortedConnMap, angle, direction))
            return cAtomParityUnknown;

        int horizontalBondType = -1;
        for (int i=0; i<mAllConnAtoms[atom]; i++) {
            if ((direction[i] & 1) == 1) { // horizontal bond
                int bondType = mBondType[mConnBond[atom][sortedConnMap[i]]];
                if (horizontalBondType != -1
                 && horizontalBondType != bondType)
                    return cAtomParityUnknown;
                horizontalBondType = bondType;
                }
            }

        // we determine the rotation for the first 3 substituents: cw/ccw
        int index = (Math.abs(direction[0] - direction[1]) == 2) ? 1 : 0;
        int dif = direction[index] - direction[index+1];
        boolean isClockwise = (Math.abs(dif) == 3) ^ (direction[index] < direction[index+1]);
        boolean is4thConnHorizontal = (mAllConnAtoms[atom] == 3 || ((direction[3] & 1) == 1));
        return (isClockwise ^ is4thConnHorizontal) ^ (horizontalBondType == cBondTypeDown) ?
                cAtomParity1 : cAtomParity2;
        }


    /**
	 * Checks whether we have two vertical non-stereo single bonds and
	 * two horizontal stereo single bonds. If these conditions
	 * are met, then an int array is returned defining the directions of all
	 * connected bonds (0:south; 1:west; 2:north; 3:east).
	 * @param atom
     * @param index map of neighbours sorted by atom index
     * @param angle bond angles sorted by neighbour atom index
     * @param direction array large enough to receive bond directions
	 * @return false if fisher projection conditions are not met
	 */
	private boolean getFisherProjectionBondDirections(int atom, int[] sortedConnMap, float[] angle, int[] direction) {
	    if (mPi[atom] != 0
	     || isAromaticAtom(atom)
	     || mConnAtoms[atom] < 3
	     || mAllConnAtoms[atom] > 4)
	        return false;

	    boolean[] isUsed = new boolean[4];
	    for (int i=0; i<mAllConnAtoms[atom]; i++) {
            float a = (float)Math.PI*5/4 - angle[i];
	        if (Math.abs(Math.PI/4 - (a % (Math.PI/2))) > FISCHER_PROJECTION_LIMIT)
	            return false;

	        direction[i] = 3 & (int)(a / (Math.PI/2));
	        if (isUsed[direction[i]])
	            return false;
	        isUsed[direction[i]] = true;

	        if ((direction[i] & 1) == 0) {  // vertical bond
	            if (mBondType[mConnBond[atom][sortedConnMap[i]]] != cBondTypeSingle)
	                return false;
	            }
	        else {
	        	if (!isStereoBond(mConnBond[atom][sortedConnMap[i]], atom))
                    return false;
	            }
	        }
	    return isUsed[0] && isUsed[2];
	    }

	
	private void setAlleneStereoBondFromParity(int atom) {
			// find preferred bond to serve as stereobond
        if (mConnAtoms[atom] != 2
         || mConnBondOrder[atom][0] != 2
         || mConnBondOrder[atom][1] != 2
         || mConnAtoms[mConnAtom[atom][0]] < 2
         || mConnAtoms[mConnAtom[atom][1]] < 2
         || mPi[mConnAtom[atom][0]] != 1
         || mPi[mConnAtom[atom][1]] != 1) {
            setAtomParity(atom, cAtomParityNone, false);
            return;
            }

        int preferredBond = -1;
		int preferredAtom = -1;
		int preferredAlleneAtom = -1;
		int oppositeAlleneAtom = -1;
		int bestScore = 0;
		for (int i=0; i<2; i++) {
			int alleneAtom = mConnAtom[atom][i];
			for (int j=0; j<mAllConnAtoms[alleneAtom]; j++) {
				int connAtom = mConnAtom[alleneAtom][j];
				if (connAtom != atom) {
					int connBond = mConnBond[alleneAtom][j];
					int score = getStereoBondScore(connBond,connAtom);
					if (bestScore < score) {
						bestScore = score;
						preferredAtom = connAtom;
						preferredBond = connBond;
						preferredAlleneAtom = alleneAtom;
						oppositeAlleneAtom = mConnAtom[atom][1-i];
						}
					}
				}
			}

        if (preferredAtom == -1)
            return;

        for (int i=0; i<2; i++)
			for (int j=0; j<mAllConnAtoms[mConnAtom[atom][i]]; j++)
				if (mConnAtom[mConnAtom[atom][i]][j] != atom)
					mBondType[mConnBond[mConnAtom[atom][i]][j]] = cBondTypeSingle;

		if (mBondAtom[1][preferredBond] != preferredAtom) {
			mBondAtom[0][preferredBond] = mBondAtom[1][preferredBond];
			mBondAtom[1][preferredBond] = preferredAtom;
			}

		int highPriorityAtom = Integer.MAX_VALUE;
		for (int i=0; i<mConnAtoms[preferredAlleneAtom]; i++) {
			int connAtom = mConnAtom[preferredAlleneAtom][i];
			if ((connAtom != atom) && (highPriorityAtom > connAtom))
				highPriorityAtom = connAtom;
			}

		int[] oppositeAtom = new int[2];
		int oppositeAtoms = 0;
		for (int i=0; i<mConnAtoms[oppositeAlleneAtom]; i++) {
			int connAtom = mConnAtom[oppositeAlleneAtom][i];
			if (connAtom != atom)
				oppositeAtom[oppositeAtoms++] = connAtom;
			}

		float alleneAngle = getBondAngle(atom, oppositeAlleneAtom);
		float angleDif = 0.0f;

		if (oppositeAtoms == 2) {
			if (oppositeAtom[0] > oppositeAtom[1]) {
				int temp = oppositeAtom[0];
				oppositeAtom[0] = oppositeAtom[1];
				oppositeAtom[1] = temp;
				}

			float hpAngleDif = getAngleDif(alleneAngle,
											getBondAngle(oppositeAlleneAtom, oppositeAtom[0]));
			float lpAngleDif = getAngleDif(alleneAngle,
											getBondAngle(oppositeAlleneAtom, oppositeAtom[1]));
			angleDif = hpAngleDif - lpAngleDif;
			}
		else {
			angleDif = getAngleDif(alleneAngle,
								   getBondAngle(oppositeAlleneAtom, oppositeAtom[0]));
			}

		if ((angleDif < 0.0)
		  ^ (getAtomParity(atom) == cAtomParity1)
		  ^ (highPriorityAtom == preferredAtom))
			mBondType[preferredBond] = cBondTypeUp;
		else
			mBondType[preferredBond] = cBondTypeDown;
		}


	/**
	 * In case bond is a BINAP kind of chiral bond with defined parity,
	 * then the preferred neighbour single bond is converted into a
	 * stereo bond to correctly reflect its defined parity.
	 * @param bond
	 */
	public void setStereoBondFromBondParity(int bond) {
        // set an optimal bond to up/down to reflect the atom parity
		if (getBondParity(bond) == Molecule.cBondParityNone
		 || getBondParity(bond) == Molecule.cBondParityUnknown
		 || !isBINAPChiralityBond(bond))
			return;

		int preferredBond = -1;
		int preferredAtom = -1;
		int preferredBINAPAtom = -1;
		int oppositeBINAPAtom = -1;
		int bestScore = 0;
		for (int i=0; i<2; i++) {
			int atom = mBondAtom[i][bond];
			for (int j=0; j<mAllConnAtoms[atom]; j++) {
				int connBond = mConnBond[atom][j];
				if (connBond != bond && getBondOrder(connBond) == 1) {
					int connAtom = mConnAtom[atom][j];
					int score = getStereoBondScore(connBond, connAtom);
					if (bestScore < score) {
						bestScore = score;
						preferredAtom = connAtom;
						preferredBond = connBond;
						preferredBINAPAtom = atom;
						oppositeBINAPAtom = mBondAtom[1-i][bond];
   						}
   					}
   				}
   			}

		if (preferredAtom == -1)
			return;

		for (int i=0; i<2; i++) {
			for (int j=0; j<mAllConnAtoms[mBondAtom[i][bond]]; j++) {
				int connBond = mConnBond[mBondAtom[i][bond]][j];
				if (connBond != bond && getBondOrder(connBond) == 1)
					mBondType[connBond] = cBondTypeSingle;
				}
			}

		if (mBondAtom[1][preferredBond] != preferredAtom) {
			mBondAtom[0][preferredBond] = mBondAtom[1][preferredBond];
			mBondAtom[1][preferredBond] = preferredAtom;
			}

		int highPriorityAtom = Integer.MAX_VALUE;
		for (int i=0; i<mConnAtoms[preferredBINAPAtom]; i++) {
			int connAtom = mConnAtom[preferredBINAPAtom][i];
			if ((mConnBond[preferredBINAPAtom][i] != bond) && (highPriorityAtom > connAtom))
				highPriorityAtom = connAtom;
       		}

		int[] oppositeAtom = new int[2];
		int oppositeAtoms = 0;
		for (int i=0; i<mConnAtoms[oppositeBINAPAtom]; i++)
			if (mConnBond[oppositeBINAPAtom][i] != bond)
				oppositeAtom[oppositeAtoms++] = mConnAtom[oppositeBINAPAtom][i];

		float binapAngle = getBondAngle(preferredBINAPAtom, oppositeBINAPAtom);
		float angleDif = 0.0f;

		if (oppositeAtoms == 2) {
			if (oppositeAtom[0] > oppositeAtom[1]) {
				int temp = oppositeAtom[0];
				oppositeAtom[0] = oppositeAtom[1];
				oppositeAtom[1] = temp;
				}

			float hpAngleDif = getAngleDif(binapAngle, getBondAngle(oppositeBINAPAtom, oppositeAtom[0]));
			float lpAngleDif = getAngleDif(binapAngle, getBondAngle(oppositeBINAPAtom, oppositeAtom[1]));
			angleDif = hpAngleDif - lpAngleDif;
			}
		else {
			angleDif = getAngleDif(binapAngle, getBondAngle(oppositeBINAPAtom, oppositeAtom[0]));
			}

		if ((angleDif < 0.0)
		  ^ (getBondParity(bond) == cBondParityZor2)
       	  ^ (highPriorityAtom == preferredAtom))
			mBondType[preferredBond] = cBondTypeUp;
       	else
       		mBondType[preferredBond] = cBondTypeDown;
		}


	protected boolean bondsAreParallel(float angle1, float angle2) {
		float angleDif = Math.abs(getAngleDif(angle1, angle2));
		return (angleDif < 0.08 || angleDif > Math.PI - 0.08);
		}


	private int preferredTHStereoBond(int atom) {
		// If we have two (anti-)parallel bonds, the we need to select
		// that one of those that is closest to the other bonds.
		float[] angle = new float[mAllConnAtoms[atom]];
		for (int i=0; i<mAllConnAtoms[atom]; i++)
			angle[i] = getBondAngle(atom, mConnAtom[atom][i]);
		for (int i=1; i<mAllConnAtoms[atom]; i++) {
			for (int j=0; j<i; j++) {
				if (bondsAreParallel(angle[i], angle[j])) {
					float angleDistanceSum1 = 0;
					float angleDistanceSum2 = 0;
					for (int k=0; k<mAllConnAtoms[atom]; k++) {
						if (k != i && k != j) {
							angleDistanceSum1 += Math.abs(Angle.difference(angle[i], angle[k]));
							angleDistanceSum2 += Math.abs(Angle.difference(angle[j], angle[k]));
							}
						}
					int bond = (angleDistanceSum1 < angleDistanceSum2) ? mConnBond[atom][i] : mConnBond[atom][j];
					if (getBondOrder(bond) == 1)
						return bond;
					}
				}
			}

		int preferredBond = -1;
		int bestScore = 0;
		for (int i=0; i<mAllConnAtoms[atom]; i++) {
			int connAtom = mConnAtom[atom][i];
			int connBond = mConnBond[atom][i];
			int score = getStereoBondScore(connBond, connAtom);
			if (bestScore < score) {
				bestScore = score;
				preferredBond = connBond;
				}
			}
		return preferredBond;
		}


	private int preferredAlleneStereoBond(int atom) {
		int preferredBond = -1;
		int bestScore = 0;
		for (int i=0; i<2; i++) {
			int alleneAtom = mConnAtom[atom][i];
			for (int j=0; j<mAllConnAtoms[alleneAtom]; j++) {
				int connAtom = mConnAtom[alleneAtom][j];
				if (connAtom != atom) {
					int connBond = mConnBond[alleneAtom][j];
					int score = getStereoBondScore(connBond, connAtom);
					if (bestScore < score) {
						bestScore = score;
						preferredBond = connBond;
						}
					}
				}
			}
		return preferredBond;
		}


	private int preferredBinapStereoBond(int bond) {
		int preferredBond = -1;
		int bestScore = 0;
		for (int i=0; i<2; i++) {
			int atom = mBondAtom[i][bond];
			for (int j=0; j<mAllConnAtoms[atom]; j++) {
				int connAtom = mConnAtom[atom][j];
				if (connAtom != mBondAtom[1-i][bond]) {
					int connBond = mConnBond[atom][j];
					int score = getStereoBondScore(connBond, connAtom);
					if (bestScore < score) {
						bestScore = score;
						preferredBond = connBond;
						}
					}
				}
			}
		return preferredBond;
		}


	/**
	 * Checks whether atom is one of the two end of an allene.
	 * @param atom
	 * @return allene center or -1
	 */
	public int findAlleneCenterAtom(int atom) {
		int center = -1;
		if (mPi[atom] == 1) {
			for (int i=0; i<mConnAtoms[atom]; i++) {
				if (mConnBondOrder[atom][i] == 2) {
					int connAtom = mConnAtom[atom][i];
					if (mConnAtoms[connAtom] == 2 && mPi[connAtom] == 2) {
						for (int j=0; j<2; j++) {
							int endAtom = mConnAtom[connAtom][j];
							if (endAtom != atom && mPi[endAtom] == 1) {
								center = connAtom;
								break;
								}
							}
						}
					break;
					}
				}
			}
		return center;
		}

	/**
	 * Checks whether atom is one of the two atoms of an axial chirality bond of BINAP type.
	 * Condition: non-aromatic single bond connecting two aromatic rings with 6 or more members
	 * that together bear at least three ortho substituents. A stereo bond indicating the
	 * chirality is not(!!!) a condition.
	 * @param atom to check, whether it is part of a bond, which has BINAP type of axial chirality
	 * @return opposite atom of axial chirality bond or -1 if axial chirality conditions are not met
	 */
	private int findBINAPOppositeAtom(int atom) {
		if (mConnAtoms[atom] == 3 && isAromaticAtom(atom) && getAtomRingSize(atom) >= 6)
			for (int i=0; i<mConnAtoms[atom]; i++)
				if (isBINAPChiralityBond(mConnBond[atom][i]))
					return mConnAtom[atom][i];
		return -1;
		}


	/**
	 * Checks whether atom is one of the two atoms of an axial chirality bond of BINAP type.
	 * Condition: non-aromatic single bond connecting two aromatic rings with 6 or more members
	 * that together bear at least three ortho substituents. A stereo bond indicating the
	 * chirality is not(!!!) a condition.
	 * @param atom to check, whether it is part of a bond, which has BINAP type of axial chirality
	 * @return axial chirality bond or -1 if axial chirality conditions are not met 
	 */
	public int findBINAPChiralityBond(int atom) {
		if (mConnAtoms[atom] == 3 && isAromaticAtom(atom) && getAtomRingSize(atom) >= 6)
			for (int i=0; i<mConnAtoms[atom]; i++)
				if (isBINAPChiralityBond(mConnBond[atom][i]))
					return mConnBond[atom][i];
		return -1;
		}


    /**
     * Checks whether this nitrogen atom is flat, because it has a double bond,
     * is member of an aromatic ring or is part of amide, an enamine or
     * in resonance with an aromatic ring. It is also checked that ortho
     * substituents don't force the amine into a non-resonance torsion.
     * State of helper arrays must be at least cHelperRings.
     * @param atom
     * @return
     */
    public boolean isFlatNitrogen(int atom) {
        if (mAtomicNo[atom] != 7)
            return false;
        if (isAromaticAtom(atom) || mPi[atom] != 0 || (mAtomQueryFeatures[atom] & cAtomQFFlatNitrogen) != 0)
            return true;
        if (mAtomCharge[atom] == 1)
            return false;
        int heteroCount = 0;
        for (int i=0; i<mConnAtoms[atom]; i++) {
        	if (mConnBondOrder[atom][i] == 1) {
	            int atomicNo = mAtomicNo[mConnAtom[atom][i]];
	        	if (atomicNo == 8 || atomicNo == 9 || atomicNo == 17)
	        		heteroCount++;
        		}
        	}
        if (heteroCount == 0) {
	        for (int i=0; i<mConnAtoms[atom]; i++) {
	            int connAtom = mConnAtom[atom][i];
	            if (mPi[connAtom] != 0) {
	            	if (isAromaticAtom(connAtom)) {
		                if (getAtomRingSize(connAtom) >= 5) {
		                    int orthoSubstituentCount = 0;
		                    for (int j=0; j<mConnAtoms[connAtom]; j++) {
		                        int ortho = mConnAtom[connAtom][j];
		                        if (ortho != atom && mConnAtoms[ortho] >= 3)
		                            orthoSubstituentCount++;
		                        }
		                    if (orthoSubstituentCount == 2
		                     || (orthoSubstituentCount == 1 && mConnAtoms[atom] == 3))
		                        continue;  // the nitrogen is rotated out of PI-plane
		                    }
		                return true;
	            		}

	            	// vinyloge amides, etc.
                    for (int j=0; j<mConnAtoms[connAtom]; j++) {
                    	if ((mConnBondOrder[connAtom][j] == 2 || isAromaticBond(mConnBond[connAtom][j]))
                    	 && isStabilizedAtom(mConnAtom[connAtom][j]))
                    		return true;
                    	}
	            	}
	        	}
        	}
        if (heteroCount < 2) {
	        for (int i=0; i<mConnAtoms[atom]; i++) {
	            int connAtom = mConnAtom[atom][i];
            	boolean isStabilized = false;
            	boolean hasCompetitor = false;
            	for (int j=0; j<mConnAtoms[connAtom]; j++) {
            		if (mConnAtom[connAtom][j] != atom) {
	            		if (mConnBondOrder[connAtom][j] != 1
	            		 && (mAtomicNo[mConnAtom[connAtom][j]] == 7
	            		  || mAtomicNo[mConnAtom[connAtom][j]] == 8
	            		  || mAtomicNo[mConnAtom[connAtom][j]] == 16))
	            			isStabilized = true;
	
	            		if (mConnBondOrder[connAtom][j] == 1
	               		 && mAtomicNo[mConnAtom[connAtom][j]] == 7)
	            			hasCompetitor = true;
            			}
            		}
            	if (isStabilized && (!hasCompetitor || heteroCount == 0))
            		return true;
                }
            }
        return false;
        }

    
    /**
	 * Checks whether bond is an axial chirality bond of the BINAP type.
	 * Condition: non-aromatic, non-small-ring (<= 7 members) single bond
	 * connecting two aromatic rings with 6 or more members each
	 * that together bear at least three ortho substituents. A stereo bond indicating the
	 * chirality is not(!!!) a condition.
	 * @param bond
	 * @return true if axial chirality conditions are met
	 */
	public boolean isBINAPChiralityBond(int bond) {
		if (mBondType[bond] != cBondTypeSingle
		 || isAromaticBond(bond)
		 || (isRingBond(bond) && getBondRingSize(bond) < 7))
			return false;

		int atom1 = mBondAtom[0][bond];
		if (!isAromaticAtom(atom1)
	     || getAtomRingSize(atom1) < 6)
			return false;

		int atom2 = mBondAtom[1][bond];
		if (!isAromaticAtom(atom2)
	     || getAtomRingSize(atom2) < 6)
			return false;

		int orthoSubstituentCount = 0;
		for (int j=0; j<mConnAtoms[atom1]; j++) {
			int connAtom = mConnAtom[atom1][j];
			if (connAtom != atom2 && mConnAtoms[connAtom] > 2)
				orthoSubstituentCount++;
			}
		for (int j=0; j<mConnAtoms[atom2]; j++) {
			int connAtom = mConnAtom[atom2][j];
			if (connAtom != atom1 && mConnAtoms[connAtom] > 2)
				orthoSubstituentCount++;
			}
		return (orthoSubstituentCount > 2);
		}


	protected boolean validateBondType(int bond, int type) {
        boolean ok = super.validateBondType(bond, type);

        if (ok && type == cBondTypeCross) {
            ensureHelperArrays(Molecule.cHelperRings);
            ok &= !isSmallRingBond(bond);
            }

        return ok;
        }

    
	public void validate() throws Exception {
        float avbl = getAverageBondLength();
        float minDistanceSquare = avbl * avbl / 16.0f;
        for (int atom1=1; atom1<mAllAtoms; atom1++) {
            for (int atom2=0; atom2<atom1; atom2++) {
                float xdif = mAtomX[atom2] - mAtomX[atom1];
                float ydif = mAtomY[atom2] - mAtomY[atom1];
                if ((xdif*xdif + ydif*ydif) < minDistanceSquare)
                    throw new Exception("The distance between two atoms is too close.");
                }
            }

        ensureHelperArrays(cHelperNeighbours);
		int allCharge = 0;
		for (int atom=0; atom<mAtoms; atom++) {
			if (getOccupiedValence(atom) > getMaxValence(atom))
				throw new Exception("atom valence exceeded");
			allCharge += mAtomCharge[atom];
			}
		if (allCharge != 0)
			throw new Exception("unbalanced atom charge");
		}

	/**
	 * Normalizes different forms of functional groups (e.g. nitro)
	 * to a preferred one. This step should precede any canonicalization.
     * @return true if the molecule was changed
	 */
	public boolean normalizeAmbiguousBonds() {
		ensureHelperArrays(cHelperNeighbours);
		boolean found = false;
		
		for (int atom=0; atom<mAtoms; atom++) {
			if (mAtomicNo[atom] == 7
			 && mAtomCharge[atom] == 0) {
				int valence = getOccupiedValence(atom);
				if (valence == 4) {
					// normalize -N(=O)-OH to -N(+)(=O)-O(-)
					for (int i=0; i<mConnAtoms[atom]; i++) {
						int connAtom = mConnAtom[atom][i];
						if (mConnBondOrder[atom][i] == 1
						 && mAtomicNo[connAtom] == 8
						 && mConnAtoms[connAtom] == 1
						 && mAtomCharge[connAtom] == 0) {
							found = true;
							mAtomCharge[atom]++;
							mAtomCharge[connAtom]--;
							break;
							}
						}
					}
				else if (valence == 5) {
					// normalize -N(=O)2 to -N(+)(=O)-O(-), -N#N to -N(+)=N(-),
					for (int i=0; i<mConnAtoms[atom]; i++) {
						int connAtom = mConnAtom[atom][i];
						int connBond = mConnBond[atom][i];
						if (mConnBondOrder[atom][i] == 2
						 && mAtomicNo[connAtom] == 8) {
							found = true;
							mAtomCharge[atom]++;
							mAtomCharge[connAtom]--;
							mBondType[connBond] = Molecule.cBondTypeSingle;
							break;
							}
						if (mConnBondOrder[atom][i] == 3
						 && mAtomicNo[connAtom] == 7) {
							found = true;
							mAtomCharge[atom]++;
							mAtomCharge[connAtom]--;
							mBondType[connBond] = Molecule.cBondTypeDouble;
							break;
							}
						}
					}
				}
			}

		// split covalent bonds between hetero atoms and one of (Li,Na,K,Mg,Ca,...)
		boolean bondDeleted = false;
		for (int bond=0; bond<mBonds; bond++) {
		    if (getBondOrder(bond) == 1) {
    		    for (int i=0; i<2; i++) {
    		        if (isElectronegative(mBondAtom[i][bond])) {
        		        int atom = mBondAtom[1-i][bond];
        		        int atomicNo = mAtomicNo[atom];
        		        if (atomicNo == 3     // Li
                         || atomicNo == 11    // Na
        		         || atomicNo == 12    // Mg
                         || atomicNo == 19    // K
                         || atomicNo == 20    // Ca
                         || atomicNo == 37    // Rb
                         || atomicNo == 38    // Sr
                         || atomicNo == 55    // Cs
                         || atomicNo == 56) { // Ba
                            mAtomCharge[atom]++;
                            mAtomCharge[mBondAtom[i][bond]]--;
                            mBondType[bond] = cBondTypeDeleted;
                            bondDeleted = true;
    		                }
    		            break;
    		            }
    		        }
		        }
		    else if (mBondType[bond] == cBondTypeMetalLigand) {
                mBondType[bond] = cBondTypeDeleted;
                bondDeleted = true;
		        }
		    }
		if (bondDeleted) {
		    compressMolTable();
		    found = true;
		    }

		if (found)
			mValidHelperArrays = cHelperNone;

		return found;
		}


	/**
	 * Canonizes charge distribution in single- and multifragment molecules.
	 * Neutralizes NH+ and appropriate number of negative hetero atoms in order to achieve a neutral molecule.
	 * This method does not change the overall charge of the molecule.
	 * @return remaining overall molecule charge
	 */
	public int canonizeCharge(boolean allowUnbalancedCharge) throws Exception {
		ensureHelperArrays(cHelperNeighbours);

		for (int bond=0; bond<mAllBonds; bond++) {
			if (getBondOrder(bond) < 3) {
				int atom1,atom2;
				if (mAtomCharge[mBondAtom[0][bond]] > 0
				 && mAtomCharge[mBondAtom[1][bond]] < 0) {
					atom1 = mBondAtom[0][bond];
					atom2 = mBondAtom[1][bond];
					}
				else if (mAtomCharge[mBondAtom[0][bond]] < 0
					  && mAtomCharge[mBondAtom[1][bond]] > 0) {
					atom1 = mBondAtom[1][bond];
					atom2 = mBondAtom[0][bond];
					}
				else
					continue;

				if (mAtomicNo[atom1] < 9)
					if (getOccupiedValence(atom1) > 3)
						continue;

				mAtomCharge[atom1] -= 1;
				mAtomCharge[atom2] += 1;
				if (getBondOrder(bond) == 1)
					mBondType[bond] = cBondTypeDouble;
				else
					mBondType[bond] = cBondTypeTriple;
				mValidHelperArrays = cHelperNone;
				}
			}

		int overallCharge = 0;
		int negativeAtomCount = 0;
		for (int atom=0; atom<mAllAtoms; atom++) {
			overallCharge += mAtomCharge[atom];
			if (mAtomCharge[atom] < 0)
			    negativeAtomCount++;
		    }

		if (!allowUnbalancedCharge && overallCharge != 0)
			throw new Exception("molecule's overall charges are not balanced");

		ensureHelperArrays(cHelperNeighbours);
		for (int atom=0; atom<mAllAtoms; atom++) {
			if (mAtomCharge[atom] > 0) {
				boolean hasNegativeNeighbour = false;
				for (int i=0; i<mConnAtoms[atom]; i++) {
					if (mAtomCharge[mConnAtom[atom][i]] < 0) {
						hasNegativeNeighbour = true;
						break;
						}
					}
				if (!hasNegativeNeighbour) {
					int chargeReduction = Math.min(getImplicitHydrogens(atom), mAtomCharge[atom]);
					if (chargeReduction != 0) {
						overallCharge -= chargeReduction;
						mAtomCharge[atom] -= chargeReduction;
						mValidHelperArrays &= cHelperNeighbours;
						}
					}
				}
			}

		if (overallCharge < 0) {
		    int[] negativeAtom = new int[negativeAtomCount];
		    negativeAtomCount = 0;
		    for (int atom=0; atom<mAllAtoms; atom++) {
				if (mAtomCharge[atom] < 0) {
					boolean hasPositiveNeighbour = false;
					for (int i=0; i<mConnAtoms[atom]; i++) {
						if (mAtomCharge[mConnAtom[atom][i]] > 0) {
							hasPositiveNeighbour = true;
							break;
							}
						}
                    if (!hasPositiveNeighbour) {
                        // Ideally priorities for negative atom protonation should
                        // be done based on atom priorities from Canonizer.
                        negativeAtom[negativeAtomCount++] = (mAtomicNo[atom] << 16)
                                                          + atom;
                        }
				    }
		        }
		    java.util.Arrays.sort(negativeAtom);
		    for (int i=negativeAtom.length-1; (overallCharge < 0) && (i>=negativeAtom.length-negativeAtomCount); i--) {
		        int atom = negativeAtom[i] & 0x0000FFFF;
		        if (isElectronegative(atom)) {
			        int chargeReduction = Math.min(-overallCharge, -mAtomCharge[atom]);
			        overallCharge += chargeReduction;
					mAtomCharge[atom] += chargeReduction;
					mValidHelperArrays &= cHelperNeighbours;
		        	}
				}
			}

		if (!allowUnbalancedCharge && overallCharge != 0)
			throw new Exception("charge canonization failed");

		return overallCharge;
		}


	/**
	 * While the Molecule class covers all primary molecule information, its derived class
	 * ExtendedMolecule handles secondary, i.e. calculated molecule information, which is cached
	 * in helper arrays and stays valid as long as the molecule's primary information is not changed.
	 * Most methods of ExtendedMolecule require some of the helper array's information. High level
	 * methods, e.g. getPath(), take care of updating an outdated cache themselves. Low level methods,
	 * e.g. isAromaticAtom(), which typically are called very often, do not check for validity
	 * nor update the helper arrays themselves. If you use low level methods, then you need to make
	 * sure that the needed helper array information is valid by this method.<br>
	 * For performance reasons there are <b>distinct levels of helper information</b>. (A higher
	 * level always includes all properties of the previous level):<br>
	 * <i>cHelperNeighbours:</i> explicit hydrogen atoms are moved to the end of the atom table and
	 * bonds leading to them are moved to the end of the bond table. This way algorithms can skip
	 * hydrogen atoms easily. For every atom directly connected atoms and bonds (with and without
	 * hydrogens) are determined. The number of pi electrons is counted.<br>
	 * <i>cHelperRings</i>: Aromatic and non-aromatic rings are detected. Atom and bond ring
	 * properties are set and a ring collection provides a total set of small rings (7 or less atoms).
	 * Atoms being in allylic/benzylic or stabilized (neighbor of a carbonyl or similar group) position
	 * are flagged as such.<br>
	 * <i>cHelperParities</i>: Atom (tetrahedral or axial) and bond (E/Z or atrop) parities are calculated
	 * from the stereo configurations.<br>
	 * <i>cHelperCIP</i>: Cahn-Ingold-Prelog stereo information for atoms and bonds.<br>
	 * <br>cHelperParities and cHelperCIP require a StereoMolecule!!!<br>
	 * @param one of cHelperNeighbours,cHelperRings,cHelperParities,cHelperCIP
     * @return true if the molecule was changed
	 */
	public void ensureHelperArrays(int required) {
		// cHelperNeighbours: mConnAtoms,mConnBonds,mPi for all atoms
		// cHelperRings: rings,aromaticity/allylic/stabilized for non-H-atoms only
		// cHelperParities: stereo parities for non-H-atoms/bonds only
		// cHelperCIP: mCanonizer, stereo parities for non-H-atoms/bonds only

		if ((required & ~mValidHelperArrays) == 0)
			return;

	    if ((mValidHelperArrays & cHelperBitNeighbours) == 0) {
			handleHydrogens();
			calculateNeighbours();

			mValidHelperArrays |= cHelperBitNeighbours;

			if (validateQueryFeatures()) {
				handleHydrogens();
				calculateNeighbours();
				}
			}

		if ((required & ~mValidHelperArrays) == 0)
			return;

	    if ((mValidHelperArrays & cHelperBitRings) == 0) {
			for (int atom=0; atom<mAtoms; atom++)
				mAtomFlags[atom] &= ~cAtomFlagsHelper2;
			for (int bond=0; bond<mBonds; bond++)
				mBondFlags[bond] &= ~cBondFlagsHelper2;

			findRings();

				// set aromaticity flags of explicitly defined delocalized bonds
			for(int bond=0; bond<mBonds; bond++) {
				if (mBondType[bond] == cBondTypeDelocalized) {
					mAtomFlags[mBondAtom[0][bond]] |= cAtomFlagAromatic;
					mAtomFlags[mBondAtom[1][bond]] |= cAtomFlagAromatic;
					mBondFlags[bond] |= cBondFlagAromatic;
                    mBondFlags[bond] |= cBondFlagDelocalized;
					}
				}

			for (int atom=0; atom<mAtoms; atom++) {	// allylic & stabilized flags
				for (int i=0; i<mConnAtoms[atom]; i++) {
					int connBond = mConnBond[atom][i];
					if (isAromaticBond(connBond))
						continue;

					int connAtom = mConnAtom[atom][i];
					for (int j=0; j<mConnAtoms[connAtom]; j++) {
						if (mConnBond[connAtom][j] == connBond)
							continue;

						if (mConnBondOrder[connAtom][j] > 1) {
						 	if (mAtomicNo[mConnAtom[connAtom][j]] == 6)
								mAtomFlags[atom] |= cAtomFlagAllylic;
							else {
								if (!isAromaticBond(mConnBond[connAtom][j])
								 && isElectronegative(mConnAtom[connAtom][j]))
									mAtomFlags[atom] |= cAtomFlagStabilized;
								}
							}
						}
					}
				}

				// propagate stabilized flags to vinylic positions
			while (true) {
				boolean found = false;
				for (int atom=0; atom<mAtoms; atom++) {
						 // for non-aromatic stabilized atoms with pi-electrons
					if (mPi[atom] > 0
					 && ((cAtomFlagStabilized | cAtomFlagAromatic)
					 		& mAtomFlags[atom]) == cAtomFlagStabilized) {
						for (int i=0; i<mConnAtoms[atom]; i++) {
							if (mConnBondOrder[atom][i] > 1) {
								int connAtom = mConnAtom[atom][i];
								int connBond = mConnBond[atom][i];
								for (int j=0; j<mConnAtoms[connAtom]; j++) {
									if (mConnBond[connAtom][j] != connBond) {
										int candidate = mConnAtom[connAtom][j];
										if ((mAtomFlags[candidate] & cAtomFlagStabilized) == 0) {
											mAtomFlags[candidate] |= cAtomFlagStabilized;
											found = true;
											}
										}
									}
								}
							}
						}
					}
				if (!found)
					break;
				}

			mValidHelperArrays |= cHelperBitRings;
			}
		}


	/**
	 * Usually explicit hydrogen atoms could be removed without changing a molecule,
	 * because the removal would just change them from explicit to implicit ones.
	 * Exceptions are hydrogen with isotop information, hydrogens not connected to a non-H atom,
	 * hydrogens carrying a custom label and hydrogens whose existence implicitly defines another
	 * atom to have an abnormal valence.<br>
	 * This method moves all simple hydrogen atoms and associated bonds to the end of the atom/bond tables.
	 * It sets mAtoms to exclude simple hydrogen atoms and mBonds to exclude bonds leading to them.
	 * mConnAtoms/mConnBonds/mConnBondOrder are neither used nor updated.
	 */
	private void handleHydrogens() {
		// find all hydrogens that are connected to a non-H atom and therefore can be implicit		
		boolean[] isSimpleHydrogen = new boolean[mAllAtoms];
		for (int bond=0; bond<mAllBonds; bond++)
			for (int i=0; i<2; i++)
				if (isSimpleHydrogen(mBondAtom[i][bond]) && !isSimpleHydrogen(mBondAtom[1-i][bond]))
					isSimpleHydrogen[mBondAtom[i][bond]] = true;

					// move all hydrogen atoms to end of atom table
		int lastNonHAtom = mAllAtoms;
		do lastNonHAtom--;
			while ((lastNonHAtom >= 0) && isSimpleHydrogen[lastNonHAtom]);

		for (int atom=0; atom<lastNonHAtom; atom++) {
			if (isSimpleHydrogen[atom]) {
				swapAtoms(atom, lastNonHAtom);

				// swap simple H flags also
				boolean temp = isSimpleHydrogen[atom];
				isSimpleHydrogen[atom] = isSimpleHydrogen[lastNonHAtom];
				isSimpleHydrogen[lastNonHAtom] = temp;

				do lastNonHAtom--;
					while (isSimpleHydrogen[lastNonHAtom]);
				}
			}
		mAtoms = lastNonHAtom + 1;

					// move all bonds to hydrogen to end of bond table
		if (mAllAtoms == mAtoms) {
			mBonds = mAllBonds;
			return;
			}

		boolean isHydrogenBond[] = new boolean[mAllBonds];
		for (int bond=0; bond<mAllBonds; bond++) {	// mark all bonds to hydrogen
			int atom1 = mBondAtom[0][bond];
			int atom2 = mBondAtom[1][bond];
			if (isSimpleHydrogen[atom1]
			 || isSimpleHydrogen[atom2])
				isHydrogenBond[bond] = true;
			}

		int lastNonHBond = mAllBonds;
		do lastNonHBond--; while ((lastNonHBond >= 0) && isHydrogenBond[lastNonHBond]);

		for (int bond=0; bond<lastNonHBond; bond++) {
			if (isHydrogenBond[bond]) {
				int tempInt = mBondAtom[0][bond];
				mBondAtom[0][bond] = mBondAtom[0][lastNonHBond];
				mBondAtom[0][lastNonHBond] = tempInt;
				tempInt = mBondAtom[1][bond];
				mBondAtom[1][bond] = mBondAtom[1][lastNonHBond];
				mBondAtom[1][lastNonHBond] = tempInt;
				tempInt = mBondType[bond];
				mBondType[bond] = mBondType[lastNonHBond];
				mBondType[lastNonHBond] = tempInt;
				isHydrogenBond[bond] = false;
				do lastNonHBond--;
					while (isHydrogenBond[lastNonHBond]);
				}
			}
		mBonds = lastNonHBond + 1;
		}

	private boolean isSimpleHydrogen(int atom) {
		return mAtomicNo[atom] == 1 && mAtomMass[atom] == 0
			&& (mAtomCustomLabel == null || mAtomCustomLabel[atom] == null);
		}

	/**
	 * Removes all plain explicit hydrogens atoms from the molecule, converting them
	 * effectively to implicit ones. If an associated bond is a stereo bond indicating
	 * a specific configuration, then another bond is converted to a stereo bond to reflect
	 * the correct stereo geometry. If the removal of a hydrogen atom would change an atom's
	 * implicit valance, the atom's abnormal valence is set accordingly.
	 */
	public void removeExplicitHydrogens() {
	    ensureHelperArrays(cHelperParities);	// to calculate stereo center parities
	    mAllAtoms = mAtoms;
	    mAllBonds = mBonds;
	    for (int atom=0; atom<mAtoms; atom++) {
	        if (mAllConnAtoms[atom] != mConnAtoms[atom]) {
	            // If we have an abnormal valence implicitly defined by explicit
	            // hydrogens, we need to explicitly define that abnormal valence!
	            int abnormalValence = getImplicitHigherValence(atom, false);

	            mAllConnAtoms[atom] = mConnAtoms[atom];

	            if (abnormalValence != -1) {
                    int newAbnormalValence = getImplicitHigherValence(atom, true);
                    if (abnormalValence != newAbnormalValence) {
                    	int explicitAbnormalValence = getAtomAbnormalValence(atom);
                    	if (explicitAbnormalValence == -1 || explicitAbnormalValence < abnormalValence)
                    		setAtomAbnormalValence(atom, abnormalValence);
                    	}
	                }
	            }
	        }
	    setStereoBondsFromParity();
		}


	private void calculateNeighbours() {
		mConnAtoms = new int[mAllAtoms];
		mAllConnAtoms = new int[mAllAtoms];
		mConnAtom = new int[mAllAtoms][];
		mConnBond = new int[mAllAtoms][];
		mConnBondOrder = new int[mAllAtoms][];
		mPi = new int[mAtoms];

        int[] connCount = new int[mAllAtoms];
        for(int bnd=0; bnd<mAllBonds; bnd++) {
            connCount[mBondAtom[0][bnd]]++;
            connCount[mBondAtom[1][bnd]]++;
            }

        for(int atom=0; atom<mAllAtoms; atom++) {
            mConnAtom[atom] = new int[connCount[atom]];
            mConnBond[atom] = new int[connCount[atom]];
            mConnBondOrder[atom] = new int[connCount[atom]];
            }

        for(int bnd=0; bnd<mAllBonds; bnd++) {
			int order = getBondOrder(bnd);
			for(int i=0; i<2; i++) {
				int atom = mBondAtom[i][bnd];
				mConnBondOrder[atom][mAllConnAtoms[atom]] = order;
				mConnAtom[atom][mAllConnAtoms[atom]] = mBondAtom[1-i][bnd];
				mConnBond[atom][mAllConnAtoms[atom]] = bnd;
				mAllConnAtoms[atom]++;
				if (bnd < mBonds)
					mConnAtoms[atom]++;
				if (atom < mAtoms) {
				    if (order > 1)
				        mPi[atom] += order + order - 2;
                    else if (mBondType[bnd] == cBondTypeDelocalized)
                        mPi[atom] = 2;
                    }
				}
			}

        for(int atom=0; atom<mAtoms; atom++)
            mPi[atom] /= 2;
        }


	private void findRings() {
        mRingSet = new RingCollection(this, RingCollection.MODE_SMALL_AND_LARGE_RINGS_AND_AROMATICITY);

        int[] atomRingBondCount = new int[mAtoms];
        for (int bond=0; bond<mBonds; bond++) {
            if (mRingSet.getBondRingSize(bond) != 0) {
                mBondFlags[bond] |= cBondFlagRing;
                atomRingBondCount[mBondAtom[0][bond]]++;
                atomRingBondCount[mBondAtom[1][bond]]++;
            	}
        	}
        for (int atom=0; atom<mAtoms; atom++) {
        	if (atomRingBondCount[atom] == 2)
                mAtomFlags[atom] |= cAtomFlags2RingBonds;
        	else if (atomRingBondCount[atom] == 3)
                mAtomFlags[atom] |= cAtomFlags3RingBonds;
        	else if (atomRingBondCount[atom] > 3)
                mAtomFlags[atom] |= cAtomFlags4RingBonds;
        	}

		for (int ringNo=0; ringNo<mRingSet.getSize(); ringNo++) {
			int ringAtom[] = mRingSet.getRingAtoms(ringNo);
			int ringBond[] = mRingSet.getRingBonds(ringNo);
			int ringAtoms = ringAtom.length;
			for (int i=0; i<ringAtoms; i++) {
				mAtomFlags[ringAtom[i]] |= cAtomFlagSmallRing;
				mBondFlags[ringBond[i]] |= cBondFlagSmallRing;

				if (mRingSet.isAromatic(ringNo)) {
					mAtomFlags[ringAtom[i]] |= cAtomFlagAromatic;
					mBondFlags[ringBond[i]] |= cBondFlagAromatic;
					}

				if (mRingSet.isDelocalized(ringNo))
					mBondFlags[ringBond[i]] |= cBondFlagDelocalized;

				if (mBondType[ringBond[i]] == cBondTypeCross)
					mBondType[ringBond[i]] = cBondTypeDouble;
				}
			}
		}


	private boolean validateQueryFeatures() {
			// returns true if hydrogens were deleted and, thus, mConnAtoms are invalid
		if (!mIsFragment)
			return false;

		for (int bond=0; bond<mBonds; bond++)
			if ((mBondQueryFeatures[bond] & cBondQFRingSize) != 0)
                mBondQueryFeatures[bond] |= cBondQFRing;

			// if an atom has no free valence then cAtomQFNoMoreNeighbours is not necessary
			// and cAtomQFMoreNeighbours is not possible
		for (int atom=0; atom<mAllAtoms; atom++)
			if (getOccupiedValence(atom) >= getMaxValence(atom))
                mAtomQueryFeatures[atom] &= ~(cAtomQFNoMoreNeighbours | cAtomQFMoreNeighbours);

			// approximate explicit hydrogens by query features
			// and remove explicit hydrogens except those with stereo bonds
		boolean deleteHydrogens = false;
		for (int atom=0; atom<mAtoms; atom++) {
			int explicitHydrogens = mAllConnAtoms[atom] - mConnAtoms[atom];
			if (!mProtectHydrogen && explicitHydrogens > 0) {
				if ((mAtomQueryFeatures[atom] & cAtomQFNoMoreNeighbours) == 0) {
					if (getFreeValence(atom) == 0)
                        mAtomQueryFeatures[atom] |= cAtomQFNoMoreNeighbours;
					else {
							// add query feature hydrogen to explicit hydrogens
						int queryFeatureHydrogens = 0;
						if ((mAtomQueryFeatures[atom] & cAtomQFNot0Hydrogen) == cAtomQFNot0Hydrogen)
							queryFeatureHydrogens++;
						if ((mAtomQueryFeatures[atom] & cAtomQFHydrogen) == (cAtomQFNot0Hydrogen | cAtomQFNot1Hydrogen))
							queryFeatureHydrogens++;

                        mAtomQueryFeatures[atom] &= ~cAtomQFHydrogen;
						if (getFreeValence(atom) <= queryFeatureHydrogens)
                            mAtomQueryFeatures[atom] |= cAtomQFNoMoreNeighbours;
						else if (queryFeatureHydrogens == 0)
                            mAtomQueryFeatures[atom] |= cAtomQFNot0Hydrogen;
						else	// queryFeatureHydrogens > 1
                            mAtomQueryFeatures[atom] |= (cAtomQFNot0Hydrogen | cAtomQFNot1Hydrogen);
						}
					}

				for (int i=mConnAtoms[atom]; i<mAllConnAtoms[atom]; i++) {
					int connBond = mConnBond[atom][i];
					if (mBondType[connBond] == cBondTypeSingle) {	// no stereo bond
						mAtomicNo[mConnAtom[atom][i]] = -1;
						mBondType[connBond] = cBondTypeDeleted;
						deleteHydrogens = true;
						}
					}
				}

			if ((mAtomQueryFeatures[atom] & cAtomQFRingSize) != 0
			 || (mAtomQueryFeatures[atom] & cAtomQFAromatic) != 0)
                mAtomQueryFeatures[atom] |= cAtomQFNotChain;

			if (mAtomCharge[atom] != 0)	// explicit charge superceeds query features
				mAtomFlags[atom] &= ~cAtomQFCharge;
			}
		if (deleteHydrogens)
			compressMolTable();

		return deleteHydrogens;
		}


	private void writeObject(ObjectOutputStream stream) throws IOException {}
	private void readObject(ObjectInputStream stream) throws IOException {}
	}
